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Theory of Biaxial Graded-lndex Optical Fiber 


1. Introduction 

1.1 Review of Previous Research 

The optical fiber has become a much studied transmission system due to its property of 
wave guidance with low loss. In recent years it has been shown that introducing anisotropies 
into the dielectric medium of the fiber produces several interesting features, such as control 
of power flow and reduction of peak attenuation near cutoff. 

Typically the analysis of wave propagation in a cylindrical dielectric waveguide such 
as an optical fiber is performed using a wave equation formulation. For the simple case 
of a step-index fiber a detailed analysis, including dispersion relations, cutoff conditions 
and mode designations, is presented by Snitzer [1], Paul and Shevgaonkar [2] present a 
similar analysis for a uniaxial step-index fiber and also perform a perturbation analysis to 
determine the modal attenuation constants. These are the only two cases for which exact 
solutions are known. 

For inhomogeneous fibers no exact solutions are known. For the case of an isotropic 
graded-index fiber several approximate analytic solution methods are available. These ap- 
proximate solutions all share the common assumption that the fiber is infinite in extent. In 
addition if the permittivity is assumed to vary slowly over the distance of one wavelength 
the wave equation formulation simplifies to an asssociated scalar wave equation. If the 
permittivity profile is parabolic the solution to the scalar wave equation can be written 
in terms of either Laguerre polynomials [3] if cylindrical coordinates are used or Hermite 
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polynomials [4] if rectangular coordinates are used. For arbitrary permittivity profiles the 
scalar wave equation can be solved using the well known WKB solution method [5], [6]. For 
parabolic permittivity profiles all three solution methods give identical results. Under the 
assumption that the fields are far from cutoff Kurtz and Streifer [7], [8] have shown that a 
solution to the full vector problem can be written in terms of either Laguerre polynomials if 
the permittivity profile is quadratic or asymptotically in terms of Bessel and Airy functions 
for arbitrary permittivity profiles which decrease slowly and monotonically. A comparison 
of the vector and scalar solutions for the quadratic permittivity profile implies the vector 
modes can be obtained by simply renumbering the scalar modes [9]. Using the renumbered 
scalar modes as a basis Hashimoto [10], [11], [12] and Bcuno [13], [14], [15] have developed 
two slighly different iterative methods which can be used to solve the full vector problem 
for an isotropic graded- index fiber. 

An alternate formulation of the problem is to write the four first-order differential 
equations for the tangential field components as a first-order matrix differential equation. 
For a step-index fiber with uniaxial core and cladding Tonning [16] has shown that the 
matrix formulation can be solved exactly in terms of Bessel functions. For isotropic graded- 
index fibers with arbitrary permittivity profiles Yeh and Lingren [17] have indirectly used 
the matrix fomulation in developing a numerical solution method based on the concept 
of stratification. Using the concept of transition matrices Tonning [18] has developed a 
numerical procedure which can be used to solve the matrix differential equation for isotropic 


graded- index fibers. 
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1.2 Outline of Proposed Research 

This thesis concerns itself with the general case of a biaxial graded-index fiber with a 
homogeneous cladding. Two methods, wave equation and matrix differential equation, of 
formulating the problem and their respective solutions will be discussed. 

For the wave equation formulation of the problem it will be shown that for the case 
of a diagonal permittivity tensor^, the longitudinal electric and magnetic fields satisfy a 
pair of coupled second-order differential equations. Also, a generalized dispersion relation 
is derived in terms of the solutions for the longitudinal electric and magnetic fields. For the 
case of a step-index fiber, either isotropic or uniaxial, these differential equations can be 
solved exactly in terms of Bessel functions. For the cases of an isotropic graded-index and a 
uniaxial graded-index fiber a solution using the Wentzel, Krammers and BriUouin (WKB) 
approximation technique will be shown. Results for some particular permittivity profiles 
will be presented. Also the WKB solutions will be compared with the vector solution found 
by Kurtz and Streifer [7]. 

For the matrix formulation it wiR4>e shown that the tangential components of the 
electric and magnetic fields satisfy a system of four first-order differential equations which 
can be conveniently written in matrix form. For the special case of meridional modes the 
system of equations splits into two systems of two equations. A general iterative technique, 
asymptotic partitioning of systems of equations, for solving systems of differential equations 
is presented. As a simple example, Bessel’s differential equation is written in matrix form 
and is solved using this asymptotic technique. Low order solutions for particular exam- 
ples of a biaxial and uniaxial graded-index fiber are presented. Finally numerical results 
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obtained using the asymptotic technique are presented for particular examples of isotropic 
and uniaxi al step-index fibers and isotropic, uniaxial and biaxial graded-index fibers. 

For purposes of comparison and verification a purely numeric solution method is also 
presented. The algorithm used by Yeh and Lindgren [17] is improved to handle the case of 
a uniaxial graded-index fiber. 
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2. Analytic Solutions 
2.1 Introduction 

Consider a circularly symmetric optical fiber with the geometry shown in figure 1 . The 
region 0 < p < a is referred to as the core and the region a < p < b as the cladding. 
The permeability of both the core and cladding is po, the permeability of free space. The 
permittivity of the cladding is eo€ c where eo is the permittivity of free space and e c is the 
relative permittivity of the cladding and is assumed to be a constant. The permittivity 
of core is €ol r where e r is the relative permittivity tensor of the core and in general is a 
function of position in the core. Also it is assumed that the radius of the cladding, 6 , is 
sufficiently large so that the fields in the cladding decay exponentially and are essentially 
equal to zero at p = 6 . This eliminates the need to impose boundary conditions at the 
air-cladding boundary. 

Consider the case where the permittivity in the core, e is given by 

(*i{p) 0 0 \ 

e(p) = €o(r{p) = «o o f 2 (p) 0 I ; (2-1) 

V 0 0 e 3 (/>) / p ^ z 

where ei(p), € 2 {p) and € 3 (p) are the relative permittivities in the />, <t> and 2 directions 
respectively. 

For time harmonic fields in a source free region, Maxwell’s equations can be written as 


V x H = ;Woe r E, 

( 2 - 2 ) 

V x E = — jfu’/zoH, 

(2-3) 

V D = 0, 

(2-4) 


V B = 0 


(2-5) 
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Figure 1 Geometry of the fiber 

where eo and po are the permittivity and permeability of free space, and u; is the angular 
frequency. The problem is to find a solution for eqs. (2-2) to (2-5) in cylindrical coordinates. 
If the z and <p dependence of the fields is given by 

p -jPz+jmt 


where (3 is the longitudinal wavenumber and m is any integer ( because the fields are periodic 
in 4> with period ), then eqs. (2-2) and (2-3) can be written in cylindrical coordinates as 


— H z + = U>€o€i E p , 

p 

(2 - 2a) 

. dH z 

-jfiHp = 

dp 

(2 - 26) 

- —Bp = Mo 

p dp p 

(2 - 2c) 

— E z + PE^ = -u>noH p , 
P 

(2 - 3a) 

dE, 

jPEp + — - = jvfoH*, 

dp 

(2 - 36) 

1 d , _ . jm _ _ 

- — (pEj,) E p = -ju>n 0 H z . 

P dp p 

(2 - 3c) 
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2.2 Wave Equation Formulation 


2.2.1 Derivation of Differential Equations 

Setting h = ZqH , where Zo = \J, JioJto is the impedance of free space, eqs. (2-2a) and (2-3b) 


can be written as the following system of equations in the unknowns E p and h#: 


m . 


ko€\E p Ph^ — h z , 
P 


P E p ktyh^ — j 


dEz 
dp ' 


( 2 - 6 ) 


where ko = u^eoMo is the wavenumber of free space. Similarly, eqs. (2-2b) and (2-3a) cam 
be written as; 


k 0 € 2 E 4> + fih p = j 


. dfx^ 


dp 1 


m 


(2-7) 


PEj, + k 0 h p E z . 

P 

Solving eqs. (2-6) and (2-7) for E p , h E$, and h p gives 


E " “ [ e 


E* = 


1 

p" 

1_ 

it 2 

ic tl 

1 

W 

K t 2 
1 


mk 0 0 dE 2 

h z -Jp 


h * ~ 4 1 f 


™P, 

jk 0 t\ 


dp y 

dE z 


dp J’ 


hp ~ P , 

*t2 L 


m/3 „ .. dh, 

E z + J ko — — 

L P d P 

-mkoi 2 c a dhz 

Ez-JP- 7 - | 

P dp J 


where 


*tn = *0 f n ~ /? 2 , n=l,2 


(2 - 8a) 
(2 - 86 ) 
(2 -8c) 
(2 - 8 d) 

(2-9) 


is the transverse wave number. Eq. (2-8) gives expressions for the transverse field compo- 
nents in terms of the longitudinal components E z and h z . 

The remaining two equations (2- 2c) and (2- 3c) can be written as 

(2- 10a) 
(2 - 106) 


+ — - — E p + jk 0 h z = 0, 
dp p P 


dh & hj, jm , 

— 1 h p — jkoi^Ez - 0. 

dp p P 
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Substitution of eqs. (2-8a,c) into (2- 10a), and eqs. (2-8b,d) into (2-10b) yields: 



m P z? | u 
jm 2 k o 


k nP 2 
m0 


k nP 

h z - ^-E' z +jk 0 h z = 0, 


k tiP 




h, - 

k nP 


jm k ()€2 TTlf3 I . c J? — n 

+ n ^ E -- W0 ‘'’° 3 - 


(2- 11a) 


(2 - 116) 


where ' = d/dp. Simplifying eq. (2-11) by collecting common derivatives of E z and h z gives 
the following 



(2 - 12a) 


(2 - 126) 


In general, eqs. (2-12a) and (2-12b) are coupled except for the case m = 0. This implies 


that the general solutions of eqs. (2-12) are of a hybrid type with both E z £ 0 and h z £ 0. 


Eqs. (2-12) can be written in a more covenient form if we make the following substitu 


tions 


and 


1 



p 

K t l _ 

p 

K t2 


1 2 ~~ 6 1 
62 — K 2 


2 — = 


4 


hi ri - « 2 


1(2 ^2 

/= 1,2 


(2-13) 

(2-14) 

(2-15) 



2,' 


#e c 


ei(«i - k 2 ) 


(2-16) 
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It is also convenient to make a change of variable from p to the normalized radius r where 
r = pja and a is the core radius. Using eqs. (2-13), (2-24) and (2-15), eqs.(2-12) can be 
rewritten as 


K + fi(r)K + A 2 gi (r)E t = p 2 {r)h' z + q 2 (r)h Zi (2 - 17a) 

K + ft(r)K + a2 9 2 {r)h t = Pl (r)E' z + qi(r)E z , (2 - 176) 


where ' = d/dr, A 2 = ( k^a ) 2 and 




«M(r) 


r €i(r)[€!(r) - k 2 ]’ 

f ,_v _ I _ f 2( r ) 

2 r e 2( r ) ~ « 2 ’ 

siM=r8hM-* ! ] 

e l( r ) 


1 - 


m 2 e 2 (r) 


92(r) = [f2(r) - k 2 
jrriK 


1 - 


A 2 € 3 (’*)[f2(r) - K 2 ]r 2 . 
m 2 


Pi( r ) = 
P2{r) = - 
?i( r ) = - 


A 2 [e a (r)-K 2 ]r 2 J’ 
«l( r ) - e 2 (r) 


r [ e a (r) - k 2 


g 2 (r) - ci(r) ' 

L f 2 (r) - k 2 
e ' 2 ( r ) 


jmK 

«i( r ) r 
jmK 
r _c 2 (r) - K 2 


?2(r) = 


JTTIK 


«i(r) 


f i( 7 ') 7 ' L f i( r ) - 


(2- 18a) 
(2 - 186) 
(2 - 18c) 

(2 - 18 d) 
(2 - 18c) 
(2-18/) 
(2-18$) 
(2 - 18 h) 


From eqs. (2-18e) through (2-18h) we can see that the differential equations become de- 
coupled for three particular cases. For so called meridional modes m is equal to zero and 
from eqs. (2-18e) through (2-18h) it can be seen that p 2 , qi and q 2 are also zero. For 
an isotropic and uniaxial step index fibers ei and e 2 are equal and constant, therefore, from 
eqs. (2-18e) through (2-18h) pi, p 2 , qi and p 2 are identically equal to zero. 
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2.2.2 Derivation of Dispersion Relation 

For the region r < 1, let the general solutions of eqs. (2-17a) and (2-17b) be given by 

E z = j4e(r), h z = Bh{r), (2 - 19) 

where A and B are constants. Using eqs. (2-8b,c) in eqs. (2-8b) and (2-8c) the tangential 
components rE ^ and rh $ can be written as 

+ (2 - 20a) 

rh t = jrBk(r) - (2 -20b) 

where e'(r) = ( d/dr)e(r ) and h'(r) = (d/dr)h(r). 

For the region r > 1, ej, e 2 , and e 3 are equal to a constant e c . Under these conditions 
eqs. (2-17a) and (2-17b) simplify to Bessel’s equations of the variable k t ar where fc t 2 = 
k^e c — jd 2 . For guided modes we require that /3 2 > k^t c and that the field be of the form 
e-' ,T as r tends to infinity, with 7 > 0. If we let -f 2 = -kf we can choose A' m (7ar), the 
modified Bessel function of the second kind, as the solution which satisfies the requirement 
of a decaying exponential. E z and h z can then be given by 

E z = C Km^ar), h z = DK^ar), (2 - 21) 

where C and D are constants and 7 2 = -k 2 t = (3 2 - k 2 0 e c . From eqs. (2-8b) and (2-8c) the 
tangential components rE $ and for r > 1 are given by 

= -^rCK^av) - ^DK'^ar) (2 - 22a) 

crp <17 

rh t = -^DK m (-,ar) + ih^CKUtar) (2 - 221.) 

v aY 07 
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where K'^ar) = dK^ar)/ d(yar). 

At r = 1 the tangential components of the electric and magnetic fields, E z , h z , E $ and 
hj must be continuous. Using eqs. (2-19) to (2-22) the boundary condition can be written 
as 


e(l) 

0 


3gfe‘(i) 


2 d) 


0 

Mi) 


-K m {^a) 

0 


ar 

jkptc TTt 

7 n 


0 

-Km{ Ifa) 


ay* 


For a non-trivial solution to eq.(2-23) the determinant 



(A\ 


/o\ 


B 


0 


C 


0 


\Dj ’ 

W 


(2-23) 


A = 


e(l) 

0 


0 

MM 


-K m (fa) 

0 


0 


a&‘'(D 


(2 - 24) 


must be identically equal to zero. 

For convenience let e = e(l), h = /*(!), e' = e'(l), h ' = h'(l), fc 2 n (l) = fc 2 n> A^ m = 
lf m (7a), and A'^ = A r ^( 7 a). By expanding the determinant and performing some algebraic 
manipulations the generalized dispersion relation is given by 


(rnpy 

1 11 

] 

1 1 I 

I 


K’ m e, e'l 

1 ^ 1 Ml 

\h ) 

(7a) 2 (Mia) 2 . 

.(7a) 2 (M2 a) 2 . 


_7 a K m (Mia) 2 e 

7 aA' m (k t 2 a ) 2 h _ 


(2 - 25) 
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2.2.3 Exact Solutions 

For an isotropic step index fiber «i, 1 2 and €3 are all equal to the constant e r . Eqs. (2-17) 
then simplify to 


e';+ ^ + 

T 


(k,af- = 0 , 


(2 - 26a) 
(2 - 266) 


where ( k t a) 2 = A 2 {e T - k 2 ) or k 2 = e T k 2 - 0 2 . E z and h z are then given by 


E z - AJm{k t ar ) , h z = B J m (ktar), 


(2 - 27) 


where A and B are constants. By substituting eq. (2-27) into the generalized dispersion 
relation given by eq. (2-25) and making use of the fact that for a step index fiber k 2 : = 


jfcjj - jfc'* gives the well known dispersion relation for a step index fiber. 

( v 


v 1 

■ 1 1 

2 


[(a 7 ) 2 + (aht) 2 J 

— 


e e KA 7°) , «C(M ' 

[07 A m (7fl) afcj Jm{kt a ) . 


1 E^(ja) + J_ JUkta) 


[fl'y Km{ 7 a ) 


(2 - 28) 


where 7 2 = /3 2 - e c ^o- 

For a uniaxial step index fiber <i = €2 # f 3 and ei and €3 are constants. Eqs. (2-17) 


simplify to 




-( k t a ) 2 - — 
L f i r 


h" + -h' z + 

r 


, , m 

(M 


= 0 

hz = 0 


(2 - 29a) 
(2 - 296) 


where (Jfe*a) 2 = A 2 ( €l - k 2 ) or k 2 = e^k 2 - (3 2 . By defining an anisotropy parameter 
p 2 = es/d, the solutions of eqs. (2-29a) and (2-29b) are given by 


E z = AJm{pk t ar ) , h z = BJ m {k t aT), 


(2 - 30) 
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where A and B are constants. By using eq. (2-30) in eq. (2-25) and making use of the fact 
that k 2 x = k\ 2 = k 2 the dispersion relation for a uniaxial step index fiber is given by [2] 


/ m3 \ 


2 r 


1 2 


+ 


fc + J'mipkta) 


\ k 0 J L( G 7) 2 ( ak t) 2 \ l a l Km(ia) ak t J m {pk t a) J 


1 K^a) + 1 J'Jkta] 


^7 ^m( 7^) ^ k t Jm{kt ®) 


(2-31) 


A representative case for both an isotropic and a uniaxial step-index fiber is presented. 
When m = 0 the solutions of the dispersion relations, either eq. (2-29) or (2-31), are either 
transverse electric, E z = 0 or transverse magnetic, h z = 0 and are designated by the notation 
TEon and TMon respectively where n = 1,2,3, .. .. When m > 0 the electric and magnetic 
fields for all solutions have components in the axial direction, i.e E z ^ 0 and h z ^ 0 and 
are therefore designated as hybrid modes. A hybrid mode is arbitrarily designated as EH 
(HE) if at some arbitrary reference point E z ( h z ) makes a larger contribution than h z ( E z ) 
to the transverse field. A less arbitrary classification scheme, which gives the same mode 
designations, based on the ratio of H z to E z at cutoff has been proposed by Snitzer[l] and 
refined by Safaai and Yip[19]. 

As an example of an isotropic step-index fiber the relative permittivities of the core and 
cladding are taken to be e r = and e c = n 2 c respectively where n T = 1.515 is the refractive 
index of the core and n c = 1.5 is the refractive index of the cladding. Figures 2 and 3 are 
plots of the normalized propagation constant, k ~ ft/ko, versus the normalized free space 
wavenumber, A = k^a for the cases m = 0 and m = 1 respectively. Two notable features 
are that the TEon and the TMon modes are essentially degenerate except close to cutoff 
and all modes except the HEn mode have a finite non-zero cutoff frequency. 


As an example of a uniaxial step-index fiber the relative permittivities in the core and 
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cladding are taken to be ei = £7 = n ?) e 3 = n \ e c = n c where rij = 1.515 is the 
refractive index of the core in the p and <f> directions, 713 = 2 is the refractive index of the 
core in the z direction and n c = 1.5. Figures 4 and 5 are plots of k versus kod for the cases 
m = 0 and m = 1. A comparison of eqs. (2-27) and (2-30) implies that the introduction of 
anisotropy into a step index fiber affects modes where E z makes the larger contribution to 
the transverse fields, i.e. TMon ^d EH mn modes. A comparison of figures 2 and 4 show 
that the TEon modes for the isotropic and uniaxial step-index fibers are identical while 
the TMom modes for the uniaxial case are displaced from the corresponding TMom for the 
isotropic case. Comparing figures 3 and 5 it can be seen that both the EH and HE modes 
for the uni axial fiber are displaced from the corresponding mode for the isotropic fiber. As 
expected the effect of the anisotropy is much more pronounced in the EH modes than in 


the HE modes. 
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2.2.4 WKB Solutions 

In order to solve eq. (2-17) for the case of an inhomogeneous fiber some simplifications 
must be made. If we assume the variation in the permittivities is very small over distances 
of one wavelength and the core is infinite in extent ( eliminates need to impose boundary 
conditions ) then the WKB method can be used. For the case of an isotropic or uniaxial 
graded-index fiber ei(r) = £2( r ) &nd eq. (2-17) simplifies into 


E" + -E z + A 2 gi{r)E z = 0, (2 - 32a) 

r 

h" + ~h' z + A 2 g2(r)h z = 0, (2-326) 

r 

where <71 and <72 are given by 

9l<r ) = <I(rj “ ** "AV' (2 -33a) 

= (2-336) 


Let E z be of the form 


E z = e jfco,/ ' (r) (2 - 34a) 

then 

E' t = jkoi>'E z , (2 - 346) 

E'J = [jk 0 r ~ kl(xk') 2 }E z , (2-34 c) 

Substituting eq. (2-34) into (2-32a) and dropping common factor of E z gives the following 
differential equation for ip{r) 

jk 0 xP" - klW) 2 + + A 2 gi = 0. 

T 


(2 - 35) 
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If gi is a slowly varying function of r then t/’( r ) can be approximated by 


V>(r) a \l>o{r) + -^i( r ); ( 2 36a ) 

V’' = V'o + (2 - 366) 

<*f = + l-V-W + p(V-i) 1 , (2 - 36c) 

tfo Kq 

V>" = < + (2-36 <f) 

«o 


Using eqs. (2-36a) through (2-36d) in eq. (2-35) gives the following equation relating 4>o and 
t/>i to the functions /i and g \ . 

+ W" - kM - 2ko^[ - m 2 

+ ^+4i+A 2 <7! =0 (2-37) 

r r 

Recalling that A 2 = (fc 0 a) 2 , if we equate hke powers of k 0 we obtain the following equations 
for t/’o and rt : 


(V’o ) 2 “ a V = °< 

jrp'o ~ WWi + l rt = 0. 

T 


Solving eq. (2-38a) gives 


tpo = i Q 



Eq. (2-38b) can be written as 




or 


M r ) = ^H r rt 0 )- 


(2 - 38a) 
(2 - 386) 


(2 - 39) 


(2-40) 
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Using eqs. (2-39) and (2-40) E t can be written as 

e ±jk 0 a J y/gi (r ) dr 

Using the same method h z can be found to be 

€ ±jk 0 a J y/gi(r)dr 

\/?[® 2 0 a ( r )] 1/4 

The mode condition for a WKB solution [20], [6] requires that 


k 0 o 


^ V9i{r)dr = ( n + n = 0, 1,2, . . . 


(2-41) 


(2 - 42) 


(2-43) 


where i = 1,2 and r* and r 2 are the turning points (zeroes) of <7;. An exact solution of 
eq. (2-43) is possible only for a small number of permittivity profiles. In general, eq. (2-43) 
must be solved numerically to determine the allowable modes. 

Consider the case where the permittivity profiles in the core are given by 


€*‘(r) = €i 1 — 2A;r a * i — 1,2,3 


(2 - 44) 


where cti is a parameter which describes the shape of the permittivity profile, 

A,= ^^ * = 1,2,3 (2-43) 

and e; is the relative permittivity at the center of the core. The value of the parameter 
Qi must be greater than or equal to one. Note that in the limit a* — > 00 the permittivity 
profile approaches the profile for a step index fiber. 

Let us consider the special case of an isotropic graded-index fiber. with a parabolic 
profile. Since ^(r) = e 2 (r) = € 3( r ) = c r( r ) eq. (2-33) reduces to 

9i{r) = 92 {r) = g{r) = e T {r) - k 2 - 


(2 - 46) 



where 
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and 



(2-47) 


A = - — - (2 - 48) 

2fr . 

For this choice of c r (r) it is possible to analytically solve eq. (2-43) to obtain the allowable 
modes. 

The turning points r a and r 2 , determined by setting $(r) = 0, are given by 



m = 0 


ri.2 



t/A* - 45 ] 


m / 0 


(2 - 49a) 
(2 - 496) 


where 


A = 


ir ~ 


K 2 


2c r A 


and 


B = 


2e r AA 2 . 


(2 - 50a) 


(2 - 506) 


When m = 0, substituting eq. (2-46) into eq. (2-43) and integrating, using rj and r 2 given 
by eq. (2-49a), results in the following mode condition 




m = 0 . 


(2-51) 


Solving eq. (2-51) for k gives 


k = 


l 

ko 



y/2^A 


(2 n + 


i) 


A 


m = 0 . 


(2 - 52) 
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Similarly, when m ^ 0 substituting eq. (2-46) into eq. (2-43) and integrating gives 


A v / 2f r A /r? + t \ \ 1 

2 V — 2 r i r 2 j* = (n+ -)* m/0. (2-53) 


Substituting for rj and r 2 from eq. (2-49b) and solving for k results in 


K = £ = Y f r - + 2 n + 1 ) m^O. 


(2-54) 


Plots of k versus Jk 0 a for the case m = 0 , 1 and 2, when n, = 1.515 and n c = 1.5 are shown 
in figures 5, 6 and 7. At the present time these WKB solutions will be designated by the 
notation WKB m „ where m and n correspond to the m and n in eqs. (2-52) and (2-53). 



Figure 6 WKB solution for an isotropic graded-indez fiber: m=0 


For the case of a uniaxial graded-indez fiber ej(r) ^ c 3 (r) and the functions 0 i(r) and 
ff 2 (r), given by eq. (2-33) are not equal. Comparing eqs. (2-33b) and (2-46) it is dear that if 
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€i(r) in eq. (2-33b) is equal to e(r) in eq. (2-46) then the solutions of the modal condition, 
eq. (2-43), when i = 2 are identical to the solutions for the isotropic case. 

Again, consider the case where the relative permittivities in the core have parabolic 
profiles. The solution of the mode condition, eq. (2-46), when i = 1 must be found by 
numerical integration. These solutions corresponding to the solutions of eq. (2-32a) for E z 
and will be designated as E mn modes. The solution of the mode condition when i — 2 are 
identical to the solutions for an isotropic graded-index fiber given by eqs. (2-52) and (2-53). 
These solutions correspond to solutions of eq. (2-32b) for h z and will be designated as H mn 
modes. Figures 9 and 10 are plots of k versus k^a for a uniaxial graded-index fiber for the 
cases m = 0 and 1 with nj = 1.515, n 3 = 2 and n c = 1.5. 

It is important to remember that the E z and h z given by eqs. (2-41) and (2-42) are 
not solutions of the complete vector problem given by eqs. (2-17) and (2-18) but are rather 
solutions of a related scalar problem given by eqs. (2-32) and (2-33). Assuming an infinite 
core, an alternative solution of the scalar problem for an isotropic parabolic-index fiber [4], 
[21] is given by 



and 

£mn = k hr ~ + 2n + l) (2 - 56) 

s o 

where m = 0,1,2,..., n = 0,1,2,..., so is the characteristic spot size of the medium 
defined as = a/ko\/2i r A, L™ is a generalized Laguerre polynomial and B mn is a modal 
constant. It can be readily seen that eqs. (2-54) and (2-56) are identical expressions for the 


propagation constant (3. 
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Assuming an infinite core and the fields are far from cutoff a vector analysis of an 
isotropic parabolic-index fiber [7], [9] gives the following for the transverse fields: 


= and Hl i] = x E? 


where 


*(*•) = *(-') I ™*, 1 


s o 


m^l 




and 


k^r - M™ + 2(n - 1)] i = 1; 


g 2 = 

r'm.n 


(2-56) 


(2 - 57) 


(2 - 58) 


kh T - 4(ro + 2n), i = 2, 

fl 0 

where m = 1, 2, 3, . . ., n = 1, 2,3, . . and the upper(lower) sign corresponds to i = 1(* = 2). 
When i = 1 it can be shown that the solutions correspond to HE modes while when i = 2 
the solutions correspond to EH modes. For the special case m — 0, meridional modes, it 
can be shown for the TEon modes that 


while for the TMon 


E n - 0 and E 6 = -* (2) 


E n = * (2) and E 6 = 0 


(2 - 59a) 


(2 - 595) 


where for both TEon and TMon modes 


4n 


0On ~ ^0 £ r 2 ' 


(2 - 60) 


Comparing the scalar solutions given by eqs. (2-55) and (2-56) with the vector solutions 

given by eqs. (2-57), (2-58), (2-59) and (2-60) it can be seen that 

for HE mn modes; 


/^vector 

Pm,n 


P%+Fn-1 fOT TE On, ™<*> and EH ™ modeS > 


(2-61) 


for m = 0, 1, 2, . . and n = 1,2, 3, . . . where /3£' n tor is given by eq. (2-58) and is given 


by eq. (2-56). 
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2.3 Matrix Differential Equation Formulation 


2.3.1 Derivation of Differential Equation 

In general the solution of eq. (2-17) is not possible except for the case of an isotropic or 
uniaxial core. A direct series solution for a more general case is not possible except for the 
case when m — 0. However, a series solution in this case still may not be possible due to 
the poles in / a (r), / 2 (r), g^r) and g 2 {r). The WKB solution of eq. (2-17) while useful for 
determining propagation constants away from cutoff is essentially the solution of a scalar 
wave equation. It also ignores the effects of electromagnetic boundary conditions and the 
effects of coupling between E z and h z 

An alternative formulation [18] is to write eqs. (2-2) and (2-3) as a set of four first 
order differential equations in terms of the tangential field components. This formulation 
preserves the vector nature of this problem and permits the use of the boundary conditions. 

Eqs. (2-2) and (2-3) can be rewritten as 


Ep 

H p 


1 

U>6o6l 


L 9 J 


1 


m _ _ _ 

— E z 4- i 

P J 


(2 - 62a) 
(2 - 626) 


and 


= juj^E^ - jPE p 

dp 

(2 - 63a) 

-^-{pEt) = jmEp - jup-opHr 

dp 

(2 - 636) 

dE / = yWoe 2 E* - jpHp 
dp 

(2 - 63c) 

— (pE<t>) = jmHp + jutotzpEz 

dp 

(2 - 63d) 
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where eqs. (2-62a) and (2-62b) represent two scalar equations and eqs. (2-63a,b,c,d) repre- 
sent four first order differential equations. Substituting eq. (2-62) into eq. (2-63), recognizing 
that Jfeo = uiy/eiojlo, Zo = y/ co/mo and k = /3/ko and making a change of variable from p to 
a normalized radius a, where a = kop = (koa,)r = A r gives 


= -j—h z + - «*)(«^), (2 - 64a) 

as se i se i 

-^-(s^) = — (m 2 - 6 i5 2 )/i, + j— (a/ty), (2 - 646) 

da aei aei 

^ = J — JS, - ~( e 2 - « 2 )(a^), (2 - 64c) 

£(-**) = --(m 2 - e 3 a 2 )£ 2 - j—(sE *). (2 - 64d) 

as s s 

Eq. (2-64) can be rewritten in matrix form as 

^ = -A(a)u, (2 -65a) 

as s 

where 

u = {E Z sE* h z sh+f (2 - 656) 


and 


A(a) = 


/ o o -i“ 

0 0 i-(m 2 - fl a 2 ) 

jm/c -jf(e 2 - K 2 ) 0 0 

V -j(m 2 — e 3 a 2 ) -jmn 0 0 / 


(2 - 65c) 


For the Bpecial case of meridional modes, m = 0, eqs. (2-64) can be separated into two 
systems each containing two equations. The first set corresponding to transverse magnetic 
modes can be written in matrix form as 


du(T - ) = Ia(™>(s)u<™) 

ds s 


(2 - 66a) 


where 



u <™> = (i. ,h t ) T 

and 

a<™>W=(.° ! 

\J€ 3« 0 / 

The Becond set corresponding to transverse electric modes can be written as 


du(TE) = Ia( te )(s)u (TE ) 

ds s 
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(2 - 666 ) 

(2 - 66c) 


where 


u< TE > = (k, 


and 


A( te )(s) = 


0 -j{(2 - K?) 


~3* 


(2 - 67a) 


(2 - 676) 


(2 - 67c) 


Eqs. (2-65), (2-66) and (2-67) can be solved by several different method depending upon 
the choice of permittivity profiles in the core. For the case of a step-index fiber, either 
isotropic or uniaxial, an analytic solution of eq. (2-65) in terms of Bessel functions [16], [18] 
is possible. This analytic solution is identical to the exact solutions given in section 2.2.3. 
For the case of an isotropic graded-index fiber an approximate method using the concept 
of transition matrices [16] cam be used. 

For the more general cases of a uniaxial or biaxial graded-index fiber these two previous 
methods are not applicable. The first method can be used in an approximate manner for an 
uniaxial graded-index fiber by ass umin g the permittivities are piecewise continuous. This 
is equivalent to the stratification technique which will be discussed in section 3. The second 
method can not be used for either a uniaxial or biaxial graded-index fiber because the 
formulation depends upon a symmetry in the matrix A(s) which is present only for the 
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isotropic case. 

What is needed is a solution method which can be used with all possible types of fibers, 
isotropic, uniaxial and biaxial, step or graded-index. One such method is the method of 
partitioning of systems of equations [22]. This method involves transforming a system 
of first order linear differential equations into a system of equations whose solutions are 
easier to find. The solution obtained using this method is valid wherever the Taylor series 
expansion for A (s) is valid. The form of the solution method presented in the following 
section is based on the expansion of A(ar) in positive powers of x in contrast to the usual 
form where the expansion is in terms of positive powers of 1/x. 

The reason for using this alternative formulation should now be readily apparent. If 
the relative permittivities are of the form given by eq. (2-44) the poles of A($) are located 
outside the fiber core in the region r > 1. The series expansion is therefore valid for the 
entire fiber core. In contrast the system obtained by writing eq. (2-17) in matrix form has 
poles in the region 0 < r < 1 whenever either ei(r) or e 2 (r) is not a constant. 
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2.3.2 Series Expansion 

Consider the following system of n linear differential equations 

— = — A(i)u(z), as x — ♦ 0 (2 - 68) 

dx *9 v ' v ' 

where u is a col umn vector, q is an integer greater than or equal to 1 and A is a N x N 
matrix given by 

OO 

A(i) = Y A n i" as x -* 0. (2 - 69) 

n=0 

We seek formal solutions of the form 

u(;r) = y(a:)x <r e A “^ 1 ^ (2 “ 70) 


where <r is a constant, 


9 + 1 

A(z) = Y 


n = 1 



m 


with A_ n = 0 for n > 0 and 


y(*) = Y ynXn asx ~ >0 

n=0 


(2-71) 


(2 - 72) 


Substituting eqs. (2-69) to (2-72) into eq. (2-68) and equating powers of i, we obtain 
equations to determine successively \ n , & and y n . 
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2.3.3 Asymptotic Partitioning of Systems of Equations 

It is possible to simplify the system of equations by transforming them into some special 
differential equations whose solutions are easier to find. Let 

u(ar) = P(r)v(x) (2 — 73) 


where u and v are column vectors and P(x) is a N x N nonsingular matrix. Using eq. (2-73), 
eq. (2-68) can be transformed into 


dv 1 x , . 
- = — B(.)v(r) 


(2 - 74) 


where 


B(*) = P(*) 


-l 


A(i)P(i) - x q 


dP{x)_ 

dx 


or 


X q ^l = A(z)P(i) - P(x)B(ar). 


dx 


(2 - 75) 


(2 - 76) 


(2 - 77) 


Choose P(x) such that B(z) has a Jordan canonical form. To do this, let 

OO 

B(z) = ^2 Bn 1 " as x — > 0, 

n=0 
oo 

P(*) = £P n z m as x — ► 0, 

n— 0 

where B n represents a Jordan canonical matrix. The left hand side of eq. (2-76) can then 
be written as 


, dPjx) 

dx 


= ^2 nPm*" +5 1 = ^2 (” ~ 9 + 1 ) P n-9+l* T 


(2-78) 


n= 1 n=q 

while the right hand side of eq. (2-76) can be written as 

OO r n 

A(x)P(x) - P(x)B(x) = ^ ^(A ( P n _/ - P/B n _/) 

n=0 M=0 


(2 - 79) 
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Using eqs. (2-78) and (2-77) eq. (2-76) can be written as 

jr (n - q + l)P„-, + i* n = £ £(A,P I - PjB n _/) 

n=q n— 0 2=0 

Equating like powers of x we obtain 

AoPo — PoBo = 0 

for z° and for x n , n > 1, 

n 

(n — q + 1 )P n — q-fl — ^ ^ (A/P n _f — P 
2=0 

where P n _ g +i = 0 for n - q + 1 < 0. Rewrite eqs. (2-81) and (2-82) as 

Bo = Pq'AoPo 


and 


n-l 


AoP n — PnBo = ( n ~ 9 + 1 )P n— q+1 — (A n _jPf P n B„_j) 

1=0 

where P 0 is chosen so that B 0 is a Jordan canonical matrix. Multiplying eq. 

the left by Pq 1 and pulling the first term out of the summation gives 
P^AoPn-Po'PnBo = 

= (n - q + l)Po 1 P»-,+l - Po'AnPo + Po'PoBn 

n-l 

- p- 1 £ ( A n _,P, - P,B n _0 


2=1 


Now define the matrices W n and F n as 


W n = Po'Pn 


and 


n — 1 


F n = P^AnPo + Po 1 ]C (A„-jP/ - P/B n _i). 


i=i 


(2-80) 
(2-81) 
(2 - 82) 

(2-83) 

(2-84) 
(2-84) from 

(2-85) 
(2 - 86 ) 


(2 - 87) 
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Eq. (2-85) can then be written as 

B 0 W n - W n B 0 = (n - q + l)W n _ g+1 + B n - F„ (2 - 88) 

If Po can be chosen so that Bo is diagonal then it can be seen from eq. (2-16) that 

(Bo)a = Af * = 1,2 , ...,N (2-89) 

where Aj is the i’th eigenvalue of A 0 - When B 0 is a diagonal matrix the expression B 0 W n - 
W n B 0 has zeroes along its main diagonal. Eq. (2-88) can be easily satisfied by setting 


(B ).. = j 1 _ ■?’ 

1 nj,J 10, i*j. 


(2 - 89a) 


and 


(W n ) it - = 0 for n > 0. 

For the particular case q = 1 eq. (2-88) reduces to 

(Bo - nI)W„ - W n B 0 = B n - F, 
where, using utf = (W n )^ and fn = (F n )jj, 


(2 - 896) 


(2-91) 


(Bo-nl) - W n B 0 

/ o 


(Ai - A 2 - n)to^ 2 (Ai - A 3 - n)w™ (Ai - A 3 - n)w 44 \ 


(A 2 - Ax - n)u > 21 0 

(A 3 - Ai - n)w 31 (A 3 - A 2 - n)w 32 0 

V (A 3 - Aj - n)u ; 41 (A 3 - A 2 - n)w 42 (A 3 - A 3 - n)w 43 


(A 2 - A 3 - n)w 23 (A 2 - A 3 - n)w 24 


(A 3 - A 3 - n)w 34 
0 


(2 - 92) 


and 


B n - F n = - 


0 

f l2 

Jn 

fl3 

Jn 

f 14 

J n 

\ 

f 21 

Jn 

0 

*23 

Jn 

f 4 

J n 


f n 31 

f 32 
Jn 

0 

f 

J n 


\/n 41 

/« 

Jn 

/« 

Jn 

0 

) 


(2 - 93) 
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Eqs. (2-91), (2-92) and (2-93) can be used to find W n . Note, if A,- - Xj - n = 0 and /A J f 0 
it may not be possible to find W n and therefore it may not be possible to find a solution. 

Using eqs. (2-86) to (2-90) the matrices B„ and P n , n = 1,2,3, ... can be found using 
an iterative procedure. After completing the desired number of iterations the matrices B(z) 
and P(z) can be approximated by series constructed from B n and P n , n = 1, 2, 3, — Since 
B(z) is a Jordan canonical matrix eq. (2-73) can be easily solved for the elements of the 
vector v(x). Then, using eq. (2-73) the solution for the vector u in the original problem 
can be found. 

As an example of this solution method let us consider Bessel’s equation 


r 2 ^-| + xj^- + (x 2 - m 2 )y = 0 
ax 1 ax 


(2-94) 


or equivalently 


d 


dx 



+ (x 2 - m 2 )y = 0. 


(2-95) 


Letting 


y-u i 


and 



(2 - 96) 


eq. (2-95) is transformed into 





(2-97) 


or equivalently 


du 

dx 


-A(i)u 

x 


(2- 98) 


where 


A(z) 




(2-99) 
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Comparing eq. (2-98) with eq. (2-68) we see immediately that q = 1. From eq. (2-99) we 
have A„ = 0 for n > 2 and 
0 1 


A-0 = ( ~2 n ) > A 1 


(2 - 100 ) 


.. „■ -■(: :)■ :) 

If m / 0 the eigenvalues of Ao are ±m and the matrices P 0 and Pq 1 can be chosen as 

p °=(i 1) “ d A) < 2 - l01 > 


bo that 


Bo = P 0 -A„F„=(™ _° m ) 


(2 - 102 ) 


From eq. (2-87) we have 


Fi = Po 1 A 1 P 0 = 0 


(2 - 103) 


and therefore B! = Wj = Pj = 0. Since Bi,Wi and Pi are all identically equal to zero, 
from eq. (2-87) 

F,-Ff>A*.- 5 L(- 1 V) 

from eq. (2-90) 


and from eqs. (2-91), (2-92) and (2-93) 

w -s(l t) 

Also, from eq. (2-86) 

P 2 = P 0 W 2 = ^-( 

4m \ — J 

v m+1 m+1 7 


(2 - 104) 


(2 - 105) 


(2 - 106) 


(2 - 107) 


Notice that both W 2 and P 2 are undefined when m = ±1 and as mentioned earlier a 
•olution may not be possible. For the moment ignore this problem with W 2 and P 2 and 
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proceed with the solution. It will be shown later that this problem can be alleviated by the 
appropriate choice of integration constants. 

Continuing the solution procedure We find F3, B3, W3 and P3 are all identically equal 
to zero. The fourth iteration of the procedure gives 


F 4 = 


b 4 = 


w 4 = 


8m 2 

1 

8m 2 

1 

8 m 2 


m+ 1 
2 

m-fl 

1 

m + 1 

0 

0 


-?) 

u 


m 

0 


(rn+l )(m+2) 


(m— l)(m— 2) 

0 


and 


P 4 = 


(m+l)(rji+2) (m— l)(m-2) 


8m 2 


(m+l)(m+2) (m-l)(m-2) 


The matrices B(i) and P(a:) can be approximated as 


(2 - 108a) 
(2 - 1086) 
(2 - 108c) 

(2 - 108d) 


B(;r) % B 0 + B 2 x 2 + B 4 z 4 
P(z) s; P 0 + P 2 * 2 + P 4 z 4 


(2 - 109a) 
(2 - 1096) 


Substituting eq. (2-109a) for B(z) in eq.(2-74) gives 


+ (B,),,*' 


Vl 


_ 1 x^_ _ 

x m 2m 8 m 2 (m + 1) 


vi 


and 


dv 2 
dx 


= \ [(B 0 ) 22 + (B 2 ) 22 z 2 + (B 4 ) 22 i 4 ]v 2 


_ i x z x7_ 

x m **" 2m 8m 2 (m - 1) 


v 2 


(2 - 110a) 


(2- 1106) 
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Solving eq. (2-110) for vj and v 2 gives 


ri(a;) = C\x m e 


_ *L fi i «- 

4m I ' tm(m+l) 


] 


and 


V2{x) — C 2X m e 


r 2 T, . r 2 
— m 4 m ^ tm(m — 1 ) 

e L 


(2- 111a) 


(2 — 1116 ) 


where C\ and C2 are integration constants and m/0. 

Notice that t>i(z) is undefined when m = - 1 and *>2(2) is undefined when m = 1. Also 
notice that the matrices P2 and P4 are undefined when m = ±1 and in addition P4 is 
undefined when m = ±2. In fact if more iterations are performed one would find that the 
matrix P 2 /c is undefined when m = ±1, ±2 , . . . , ±k. It appears the eventual solution for 
Ui(a*) and 1x2(2*) will always be undefined w T hen m is an integer in the interval [-k . . . k] 
where 2 k is the number of iterations performed. This problem can be easily overcome by 
setting C 2 = 0 when m > 0 and C\ = 0 when m < 0. However, the two solutions which are 
obtained are not independent solutions when m is an integer. A careful inspection of the 
expressions for tJ 1 (x), t> 2 (z), P2 and P 4 show that 


and 


m*) 


= t'i(x) 


(2 - 112a) 


(P*)« 


= (ftJa 


(2 - 1126) 


(P *)* 


(P*)a 


(2- 112c) 


where i - 1,2 and (P n ).j is an element of P n . Therefore, it is only necessary to consider 


the solutions for m > 0. 
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We can then write 


( Ul A _ ( (^ > °)n + (^Mn* 2 + (^Mn* 4 \ ^ 
{uJ-{ (P 0 ) 21 + (P 2 ) 21 * 2 + (P 4 ) 21 * 4 ) 1 


= Cjx T 


^ 4m(m+l ) 1 2m(m+2)] \ | 1 +» m (m+l)] 


m ~ 77 mTT) [l + — 


). 


(2-113) 


4(m+l) [ A ^ 2m(m+2)J 

The constant C\ can be determined by examining the solution of u 2 at x = 0 when m = 1. 


When m — 1, t* 2 (x) = x d,Ji/dx or 


dJi(z) _ u 2 (ap _ 
dx x 1 


1 - 


8 


K)] 


,-*M) 


but at x = 0, dJ a /dx = 1/2 and t* 2 /x = Cj therefore Ci = 1/2. The solution of eq. (2-94) 
for m > 0 can be written as 


y(* 


>=H 


i + 


i + 


4m(m + 1) [ 2 m(m -f 2) 


1 - — 

4m 

F 


1 + 


•m(m+l ) 


(2 - 114) 


Now consider the solution of eq, (2-94) when m = 0. From eqs. (2-98) and (2-99) we 
have <7=1, A n = 0 for n > 2, 


Aq = 






(2 - 115) 


and Ao has two eigenvalues equal to zero. Since Ao is already in Jordan canonical form let 
P 0 = P" 1 = I where I is the identity matrix, so that 


Bo = P 0 1 AqPo = 



(2 - 116) 


ig also in Jordan canonical form. Performing four iterations of the solution procedure results 
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in the following matrices: 

Fi = Bj = Wj = Pi = 0 



F 3 = B 3 = W 3 = P a = 0 - 



The matrix B(x) can be written as 

B(x) » Bo + B a x 4 = _^) 

Eq. (2-74) can then be written as 


dv i 1 

1 X = z 1 ’ 2 


and 


dv 2 
dx 



Solving eq. (2-118) first for vi and then for v\ results in 


v 2 (x) = C ie -* 4 ' 16 


m(x) = J^-e~ xi/l6 dx 

= C 2 + Cie" 1 "/ 16 ln(z) + C a J ^-e~ xt ' ie hi{x)dx 


(2 - 117) 


(2 - 117) 


(2 - 118a) 


(2- 1186) 


(2 - 119a) 


(2-1196) 


and 
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where C\ and Ci are integration constants. In order for vi(z) to be finite at x = 0 ( assumes 
the desired solution is Jo not Ko ) the constant Ci must be identically equal to zero, .or 
vi(x) = C 2 and V 2 {x) = 0. The solutions of eq. (2-98) can then be written as 


u(x) = P(x)v(x) 



When m = 0 the solution of eq. (2-94) can be written as 

»<*> = Ci (' • 7 + s) 

which is simply the tnmcated series expansion for Jq(x). 


(2 - 120 ) 


(2 - 121 ) 


2.3.4 Solutions for Transverse Modes 

Assuming the individual elements of the matrices A^ TE ^(s) and A^™^(s) can be expanded 
as Taylor series, a completely general solution to eqs. (2-66) and (2-67) can be found in 
terms of the coefficients from the series expansions. After two or three iterations of the 
solution procedure the resulting matrices become cumbersome and further iterations are 
tedious. If the form of the permittivity profiles is known in advance the iteration procedure 
can often be made more manageable. 

Let us ass um e the permittivity profiles are of the form given by eq. (2-44). In particular 
choose all the profiles to have a parabolic shape i.e. a; = 2 i = 1,2,3. It should be noted 
that since e a (r) and e 2 (r) must be equal at r = 0, it is necessary for 6j(0) = € 2 ( 0 ) and 


Ai = A 2 . Since ai = 02 = 2 this choice for the permittivity profiles does not strictly 
contain an example of a biaxial graded-index fiber. If however, in the final result either Ai 
or A 2 but not both is set to zero then the resulting solution is a valid example of a biaxial 
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graded-index fiber where either €i(r) has a parabolic profile and € 2 (r) is a constant or €i(r) 
is a constant and €2( r ) has a parabolic profile. 

If all three permittivities have parabolic profiles then 

€i(s) = €i(l - 2A °s 2 ) i = 1,2, 3 (2 - 122) 


where s = k 0 p = k 0 ar and A° = A i/(k 0 a) 2 . Substituting eq. (2-122) into the expression 
for given by eq. (2-66c) and expanding A^™) as a series results in Aj™i = 0 for 

n = 1,2,3,... and 


a (tm)_/0 £ k m 
0 0 


,(TM) _ 


0 -£(2A?)k 2 


(2 - 123) 


(TM) _ 


(TM) 
l 2 n 


0 

0 -£(2A°) 2 k 2 

-;2e 3 A° 0 ) 

_(0 -£(2A°)"k 2 \ 

0 0 ) 

for n = 1,2,3,.. . where JfcjL = €, - k 2 . Similarly the expansion of A< TE ) gives A = 0 


for n = 1 and n > 2, 


, (TE) _ 


o A,= (-°. *T°) < 2 - 124 > 


where Jfejy 2 = e 2 - * 2 . Since the two eigenvalues of both Aq™^ and Aq TE ^ are identically 
equal to zero the matrix Po in both cases must be chosen so that Bo is Jordan canonical 
matrix. Since Aq™^ and Aq TE ^ are of the form 

A„=(" l) (2-125) 


if Po is chosen as 


■(; :) 


(2 - 126) 
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then 


B " = Po- , A»P„=(o o) 


(2-127) 


where eq. (2-127) is valid for both Aq™^ and Aq TE \ 


Since Bo is not a diagonal matrix the choices made in eq. (2-89b) for the elements of 
W n will not in general satisfy eq. (2-88). For the particular case q = 1 if the elements of 
B n are chosen using eq. (2-89a) the elements of W n must be chosen so that eq. (2-91) is 
satisfied. For a 2 x 2 matrix eq. (2-91) can be written as 


- nw]^ w 22 - nw^ 2 - «' n 


— nw 


21 


-nw. 


22 


- *)“(-]? -f) (2 - 128) 


where w% = (W „)„ and $ = (F„)ij. Solving eq. (2-128) for the elements of W n gives 

/ 21 

J n 


W*J = 


n 

21 f 21 


w n = <L = f, 


22 

W = — 


n n 4 

u ; 21 -/ 21 

_ Jn 


(2 - 129) 


n n* 


.12 _ j /» _ r?Jl - 2ff 

U>n n n 3 


Eqs. (2-66) and (2-67) can now be solved iteratively using eqs. (2-86), (2-87), (2-89a) and 
(2-129). 

After performing N iterations the matrix B(s) for either eq. (2-66) or (2-67) can be 
written in the form 

*!(.)) ( 2 - i3o °) 


where 


N 




n=l 


i = 1,2. 


(2 - 1306) 
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Note that the summation in eq. (2-130b) starts at n = 1 because (B n )^ = 0 i = 1,2. 
Substituting eq. (2-130) into eq. (2-74) gives 

= -[5n(«)vi + V 2 ] (2 - 131a) 

as s 

^ = ~B 22 (s)v2 • (2-1316) 

as s 

Solving eq. (2-131) first for t> 2 (s) then for vi(s) gives 


1 . 2 ( 5 ) = C 2 e Xi{,) (2 - 132a) 

t>!(s) = + J ds (2 - 1326) 

where 

A i(s) = j-Bn($)d$ i— 1,2 


N 




(2 - 132c) 


n=l 

In order for ^i(s) to be finite at $ = 0 it is necessary to choose C 2 to be identically equal 
to zero in eq. (2-132). The solutions for i>i(s) and t> 2 (s) can then be written as 


t>i(s) = Cie A,(,) 

V 2 (5) = 0 . 

The solution to the original problem now is written as 

U(5) = P(5)v(5) 



(2 - 133a) 
(2 - 1336) 


(2 - 134a) 


where 
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N 


P.M = E < p »)^" • (2 “ 134i) 

n=0 

/ T'Vf 

If four iterations are performed for the transverse magnetic case using the A„ *s given 
in eq. (2-123) the solutions for E z and sh ^ can be written as 

Ez = f3 k2ffl - 2 


^3 


€i 4 


shj) = —Ci ■ 


„2 

ei 2 


s" + 


+ + (2 -135a) 

[Vfi/ 64 ^i 8 J I 

+ f* A oSi] j4 l e a(AX)4 . ( 2 - 1356) 

Vej/ 16 €i 2 J J 


Similarly for the transverse electric case h z and 5^ are found to be 


h z — — i^ 2 c i 
rjb 2 


4 


1- ^5 2 + 


b 4 

k N 

64 




= —Ci 


JV2 2 *A^2 4 

5 ^ i 


c «» a St 


(2 - 136a) 

- (2 - 1366) 

2 16 1 

An important question to ask at this time is whether or not these asymptotic solutions 
correspond to any known solutions, preferably an exact solution. The only comparison 
which can be made with an exact solution is for the case of either an isotropic or a uniaxial 
step-index fiber. The asymptotic solutions for the transverse modes in a step-index fiber 
are obtained by setting = 0, i = 1,2,3 and €2 = ei in eqs. (2-135) and (2-136). For the 
transverse electric case the asymptotic solutions for h z and sE $ are given by 

* -%** + &**) 


and for the transverse magnetic case E t and sh^ are given by 


E z = ^-k 2 m Ci 


(2 - 137a) 
(2 - 1376) 


[1 C3 **v + (^ 

3\ 2 *tl a 4 

61 4 v< 

:i / 64 


(2 - 138a) 
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sh^ — E 


L f3 (- 

j 4 

. f! 2 v< 

7/ 16 


(2 - 1386) 


where Jbjyj = ei — /c 2 for both eq. (2-137) and (2-138). Prom eqs. (2-8b), (2-8c) and (2-30) 
the exact solutions for the transverse electric case are given by 


„ . $ dh z . s . x 


and for the transverse magnetic case the exact solution is given by 


H z 



u ■ se ' 
S — •? t 2 
*JVl 


dE z 

ds 


D 7 / 

r — 

ks\ 



(2 -139a) 
(2 - 1396) 


(2 - 140a) 
(2 - 1406) 


If Ci = jB/k 2 Nl in eq. (2-137) and C x = -je l Bjk 2 Nl in eq. (2-138) then it is clear that 
eq. (2-137) and (2-138) are simply truncated series expansions for eqs. (2-139) and (2-140). 

The allowable modes can be determined by solving the generalized dispersion relation 
given by eq. (2-25) For transverse modes the generalized dispersion relation separates into 
the following two equations 


1 K'Jna) 1 V = Q 

7 a -Km (7 a ) (k t2 a) 2 h 

KnM £l j _ Q 

7 a K m (~ta) (k t2 a) 2 e 


(2 - 141a) 
(2 - 1416) 


where e is the solution for E z in the core evaluated at r = 1 (s = koa), e' = de/dr , h is the 


solution for h z in the core evaluated at r — 1, h' — dh/dr, kti and kt 2 arc the transverse 
wavenumbers, eq. (2-9), evaluated at r = 1. Eq. (2-141a) is the dispersion relation for 
transverse electric modes and eq. (2-141b) is the dispersion relation for transverse magnetic 
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modes. For transverse electric modes the ratio h! fh is given by 


h!_ 

h 


dhj 

dr 

h z 


kod^f 

___ ’^t2 { s ) sE<f> 

K 

~ 3 U2 h 

"-O UZ 

8=koa 


$=koa 


= -3 


k? 2 (k 0 a) P 21 (fc 0 a) 


(2 - 142) 


kl P\\{k 0 a) 

where k\ 2 is given by eq. (2-9). For transverse magnetic modes the ratio t' j e is given by 


e E z 


koa^t 

F 2 j(s) sh 

,,r e ■ 

" fc^ei(s) E z 

8 ~kq a 


»=fcoa 


_ k^(k 0 a) P2i{k 0 a) 
k%€i(k 0 a) Pn{k 0 a) 

where k ^ is given by eq. (2-9). 


(2 - 143) 


2.3.5 Solution for Hybrid Modes 

As was the case for the solution of eqs. (2-66) and (2-67) for the transverse modes it is 
possible to generate a general solution for eq. (2-65) in terms of the elements of the coefficient 
matrices from the series expansion of A(s). In practice, however, it is not desirable to 
generate such a solution. Instead for mathematical convenience consider the solutions of 
eq. (2-65) for some particular permittivity profiles. 

The two example which will be discussed are a biaxial graded-index fiber where e 2 (s) 
is a constant and a uniaxial graded-index fiber. In both cases €i(s) and 63(5) ax e chosen 
to have parabolic profiles. These two examples contain as special cases the solutions for a 
step-index fiber, either isotropic or uniaxial, and an isotropic graded-index fiber. 

For the example of a biaxial graded-index fiber, if ei(s) and €3(5) are given by eq. (2-122) 
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the matrix A(s) can be expanded in terms of the following matrices 


/ 


An = 


0 

0 


k? ^ 
_ jtaSi jtm. \ 


A,= 


jmn 

\ -jm 2 

0 

0 j 

0 


~3 «i 3 e , 

imi 

0 


■mK 
J <1 

0 


A 4 = 


\-;M 2A°) 


/ 0 

0 
0 

Vy«3 o 

o 

0 
0 
0 


0 
0 

—jmK 0 0 / 

0 0 


(2 - 144a) 


0 


0 


(2 - 1446) 


/ 


-j^(2A°) 2 -j^-(2A°) 2 \ 
^(2A?) 2 J^(2A?) 2 


0 

0 


0 

0 


(2 - 144c) 


/ 


0\n 


A2n = 


(2A?) 


/0 0 
0 0 
0 0 
\0 0 


-jmK —jK l \ 

,2 


jm 

0 

0 


JTTIK 

0 

0 


n = 3,4,5, . . 


(2 - 144d) 


and A 2n -i = 0 for n = 1,2,3, . . . where k 2 Nl = ei - k 2 . Note, e 2 does not appear in eq. (2- 
144) since « 2 (0) = ei(0) = For the example of a uniaxial graded-index fiber where £a(s) 


and f 3 (s) are again given by eq. (2-122) the matrix A(s) can be expanded as a series using 

/ 0 0 


An = 


* TTIK - 

-3 IT J 


\ 


0 0 
jmn 
\ - j™ 2 


A 2 = 


/ 0 

0 

0 

3 


o 

0 

* i(2A?) 


-jmn 

-J“( 2A?) 

.r 

j 


jm 

ci 

o 

o 


* rnK 
3 — 

0 


(2 - 145a) 


0 ) 


~3\ f(2A?)\ 


A 4 = 


0 
0 

0 0 

V-je 3 (2A«) 0 


£ (2A?) - l] 

0 
0 

0 -j*f(2A°) 2 -j£(2A°) 2 ^ 

0 *^(2A°) 2 j^( 2A?) 2 


J^(2A?) 

0 

0 


(2 - 1456) 


0 

0 


(2 - 145c) 
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^2n = 



(0 

0 

—jmK 

-j« 2 \ 

r 

0 

0 

jm 2 

jmK 


0 

0 

0 

0 


Vo 

0 

0 

0 / 

1 , 2 , 3 ,. 


where : 

= ~ 


n = 3,4,5, 


(2 - 145d) 


it can be readily seen that the series expansion for both cases are identical except for the 
presence of an additional element in A 2 for the uniaxial graded-index case. The eigenvalues 
of A 0 for both cases are ±m. Since the A 0 ’s have repeated eigenvalues, in general, the choice 
for the matrix P 0 should at best cause B 0 to be a Jordan canonical matrix. Note, this is 
the only restriction which the solution method places on the form of P 0 . Any P 0 which 
causes B 0 to be a Jordan canonical matrix can be expected to result in a valid solution. 


Since it is posssible for several different choices of P 0 to satisfy this condition, conceivably 
there may exist several possible mathematical solutions to the problem. 


Since the solution for a step-index fiber exists as a special case of the solution for a 
graded-index fiber it is reasonable to choose P 0 based on the knowledge of the exact solution 
for a step-index fiber. For the case of a uniaxial step-index fiber the exact solutions for E z 
and h z as given by eq. (2-30) suggest that Po should be chosen so that the resulting P( 5 ) 
yields E z = P u v j + P 13 v 3 and h z = P 37 v 2 + P 34 v 4 as solutions. If P 0 is chosen as 



JU2 

K N\ 

0 

U 2 
K N 1 

0 ^ 


rriK 

jm 

mn 

— jm 


0 

1 

0 

h 2 

k n\ 


-jme x 

mK 

jmti 

mK / 


(2 - 146) 


this additional requirement is at the least satisfied for the lowest order solution where 


P(s) = Po- 
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Using eq. (2-102) B 0 is given by 

B 0 = 


( m 

0 

0 

° \ 

0 

m 

0 

0 

0 

0 

—m 

0 

o 

0 

0 

— m / 


(2 - 147) 


This is the more convenient form for Bo since eq. (2-88) can be easily satisfied by choosing 
B n and W n according to eq. (2-89) and hence W n can be found using eqs. (2-92) and 
(2-93). If Bo was not a diagonal matrix then W n would have to be found from a more 
complicated expression similiar to eq. (2-128). 

Since Bo is a diagonal matrix, choosing B n according to eq. (2-89a) makes B(s) a 


diagonal matrix. Therefore, in general, B(s) can be written as 

(B n {s) 0 0 0 \ 

B (*) = 


0 ^22( 5 ) o o 

0 0 I?33(s) 0 

V 0 0 0 B AA (s)f 


(2 - 148a) 


where 


N 


Bu{s) = ±m + ]T(B n ) ii5 n * = 1,2,3,4 , (2-1486) 

n — 1 

N is the number of iterations and the upper(lower) sign corresponds to i — l,2(t = 3,4). 
Using eq. (2-74) the differential equation for v(s) can be simply written as 


= -Bu(s)ds i = 1,2, 3, 4. 
as s 


(2 - 149) 


The solution of eq. (2-149) for t>, is then given by 


Vi (s) = CiS ±m e X ' {e) i = 1,2, 3, 4 (2 - 150a) 

where Ci i = 1, 2, 3, 4 is a constant, 

Ms) = j - s [Bu{s)^m}ds i = 1,2, 3, 4 

" ,n 

n=l 


*■ = 1 , 2 , 3, 4 


(2 - 1506) 
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and the upper(lower) sign corresponds to i = 1,2 (t = 3,4). Since 11 ( 5 ) must be finite at 
5 = 0 v(s) must also be finite at s = 0 and it is therefore necessary to set C$ = C 4 = 0. 
The solution to eq. (2-65) can then be written as 


Ez \ 


/Pu{») Pni s ) \ 

sE$ 

- s m 

P2l{s) ^22(5) 

K 

— 0 

P 31 {s) ^32(5) 

V sh.4 } 


\p*i(s) pm) 


/Cie A >W\ 

\C 2 e x ’(‘)J 


(2- 151a) 


where 


N 


P ij (s) = 'jT( P "h sT> * = 1,2,3,4; j = 1,2, 3, 4. (2-1516) 

n=0 

For the example of a biaxial graded-index fiber the following expressions for A,(.s) and 
Pij{s) are obtained after two iterations 


M*) = 

Pl \{ s ) = («1 - K 2 ) + 




4m(m + 1) 
^ . . .m/c /m + 2 \ . A - 

P^s) = - 3 - r (— T )(2A?)-- 


(f! - K 2 ) 2 - m 2 /c 2 (2Ai)| s 2 


(2 - 152a) 
(2 - 1526) 

(2 - 152c) 
(2 - 152d) 


P 21 (& ) = mK + 


4(m 


K f f 3 , 2 


+ [(ci - /c 2 ) -f (m + l)ti )s 2 (2 — 152e) 


= - K,) 


_ , . . mrt\K /rt a0 v 2 

J , 3i(^)--J 4(m + 1 j(2A 1 ) 5 

-P32(«) = («1 - « 2 ) + 


+ m2(2A | - f(m + 1)/C 2 - (d - k 2 )} Is 2 (2 - 152/) 
— /C L J J 

(2 - 152<?) 


4m(m + 1) 


[(«i - k 2 ) 2 - m 2 e a (2A?) . 


(2- 1526) 
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P 42 (S) = ^+ 4(m ' t +1) [( € l - * 2 ) 1 ( 2A ?)] 52 ' 


(2 - 152z) 
(2 - 152j) 


This should not be considered an accurate solution for u(s) since the term A3 does not 
appear anywhere in eq. (2-152). This solution is identical to the solution obtained after two 
iterations for a biaxial graded-index fiber where (s) has a parabolic profile and 62 (5) and 
63(5) are constant. Since A^ only appears in the matrix A 4 at least four iterations must be 
performed in order to obtain the effects of a non-constant €3. A solution for a uniaxial or 
a step index-fiber can be obtained from eq. (2-152) by setting Aj( and A3) equal to zero. 

For the example of a uniaxial graded-index fiber two iterations produce the following 
expression for A t (,s) and pm 


A 1 ( i ) = --Lli( fl -K 2 )5 2 (2 -153a) 

4m €1 

A ; (.) = - «V (2 - 1536) 


Pn (.) = (<1 - «’) + 4m( ^ +1) [|(«. - «’)' - 2mV(2A”)]. J (2 - 153c) 


**<•> = (2 " 16M) 
/■„(.) = m« + — ^ [|(<, - *’) + 2m»(2A»)] (2 - 153e) 

i>22« = jm - 4 ,J + t) [(«i - «*) - 2m I (2A°)j « ! (2 - 153/) 

™ = (2 - 153S) 


■P32(j) = ( f l ~ r2 ) + 


(2 - 153 h) 



52 


Pai(s) = - jme i + ~ k 2 )* 2 (2 - 153t) 

4 (m + 1) 

p 42 ( 5 ) = + - ~ (e l - k 2 )s 2 (2 - 153;) 

4(m + 1) 


As was the case for the example of a biaxial fiber, the expressions in eq. (2-153) are not an 
accurate solution since the term A® does not appear in any equation. Again, in order to 
see the effects of €3(5) it is necessary to perform at least four iterations. 

For the solutions of eq. (2-65) the generalized dispersion relation, given by eq. (2-25) 
can not be used. For the hybrid modes it is possible to derive from eq. (2-64) expressions 
for e'/e and h* /h' similar to eqs. (2-142) and (2-142). However, in general e'/e and h! jh 
are functions of the unknown constants C-\ and C 2 which appear in the general solution 
for n(s). For special cases, such as a step-index fiber, where E z and h z are decoupled it 
may be possible to set either C\ or C2 equal to zero without losing a a complete solution. 
The generalized dispersion relation can only be used if either C\ or C 2 can be set to zero. 
Instead, using the solutions for eq. (2-65) a new dispersion relation must be derived by 
enforcing the electromagnetic boundary conditions at the core-cladding interface. 

Using eq. (2-21), (2-22) and (2-151) the boundary conditions at s = k 0 a is satisfied 
provided 


/Pll Pi 2 -Km 0 \ 


Ci{k 0 a) m e Xl \ 


(°\ 

P21 P22 ^K m j^K' m 


C 2 {koa) m e Xi 


0 

Pz\ Pz2 0 —Km 


C 


0 

\P« p,2 WK m ) 


\D ) 


w 


(2 - 154) 


where Pij = Pij(k 0 a), A, = A<(feoa), K m = K m {koa~f N ), K' m = K^koa^s) and 7^ = 
k 2 - e c . A non-trivial solution of eq. (2-154) exists whenever the determinant is equal to 


rero. An explicit equation for the determinant is not provided since it cannot be expressed 
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in a convenient form as was the case for the generalized dispersion relation derived in section 
2.2.2 . Instead the zeroes of the determinant are found directly from eq. (2-154). 

2.3,6 Numerical Results 

As was previously stated the solutions for the biaxial graded-index fiber and the uniaxial 
graded-index fibers given by eq. (2-152) and (2-153) do not include the effects of a non- 
constant e 3 (s). Obtaining a more accurate solution requires performing more than four 
iterations. For the example problem shown in section 2.3.3 performing more than four 
iterations is feasible since A n = 0 when ti > 2 and A(s) is a 2 x 2 matrix. In contrast, 
performing more than two iterations in order to solve eq. (2-65) is much more difficult since 
A (s) is a 4 x 4 matrix and in general A n ^ 0 for n > 2 when ei(s) is not constant. 

Instead of deriving algebraic equations for the elements of F n , B n , W n and P„ the 
values of these matrices can be determined numerically if the values of m, k and k 0 a are 
known in advance. There are two difficulties with this appoach. First, numerical errors can 
develop since the accuracy of the matrices obtained in the t’th iteration depends upon the 
accuracy of the matrices obtained in the previous i - 1 iterations. The second and more 
important problem comes from method in which the matrices B„ and W n are chosen. As 
previously mentioned, it may not be possible to find W n whenever A; — A j — n = 0 where A, 
and Xj are eigenvalues of A 0 . In the analytic solution one could simply ignore the problem 
during the iteration process and then at the end of the process throw out the unbounded 
solutions with an appropriate choice of constants. The ability to do this appears to depend 
upon the form of A(a) and the ordering of the eigenvalues of A 0 in B 0 Luckily, due to 
the form of A(a) setting the third and fourth columns of W n equal to zero is equivalent to 



54 


setting Cz and C 4 equal to zero as was done in the solution of v($) given by eq. (2-150). 

Figures 11, 12 and 13 Eire plots of k versus koa when m = 1 for the examples of 
the five possible types of fiber cores previously discussed. In each example e\ = n\ and 
c c = t? c where nj = 1.515 and n c — 1.5. For the uniaxial fibers, step-index and graded- 
index, 63 = n\ where n 3 = 2.0. For the isotropic and uniaxial graded-index fibers all the 
permittivity profiles are parabolic. For the biaxial graded-index fiber ei and e 3 are parabolic 
while €2 is a constant. Due to the form of the solutions when m/0 only the mode with 
the lowest cutoff frequency can be found for a given value of m, which for m = 1 is the 
HEn mode. Based upon the extremely poor agreement between the asymptotic and exact 
solution method for a step-index fiber no results are given for the solutions of eqs. (2-66) 
and (2-67). For the case when m = 0 the asymptotic solution for a step index fiber is 
equivalent to finding a series solution for u(s). The poor agreement can then be attributed 
to using too few terms in the series to approximate the solution and can also be due to 
numerical errors from the evalution of the series. 

In both figures 11 and 12 the asymptotic solutions for the isotropic and uniaxial step- 
index fibers are in good agreement with exact results. For the isotropic step-index fiber the 
HEn mode for the asymptotic and exact solutions are almost identical when koa < 10. For 
Jboa > 10 the asymptotic solution begins to diverge from the exact solution but for koa > 20 
the distance between the two curves is approximately constant. For the uniaxial step-index 
fiber the asymptotic solution begins to diverge from the exact solution around koa = 8 but 
the separation distance is reasonably constant for k 0 a > 20. 


Figure 13 is a plot of k versus koa for an isotropic, a uniaxial and a biaxial graded-index 
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fiber. As was the case for the step-index fiber the HEn mode for the uniaxial graded- 
index fiber is slightly displaced from the HE n mode for the isotropic case. However, the 
displacement for the uniaxial graded-index fiber is not as large as the displacement for the 
uniaxial step-index fiber. A comparison of figure 13 with either figures 3 and 5 or figures 
11 and 12 shows that for a given value of k 0 a the value of k for the HEn mode in either 
an isotropic or uniaxial graded-index fiber is less than the value for the corresponding step 
index fiber. Also for the HEn mode k as a function of koa for an isotropic or uniaxial 
graded-index fiber increases less rapidly than in a step-index fiber. This indicates the pulse 
delay which is proportional to d.K/d{koa ) is smaller for the HEn mode in either type of 
graded-index fiber than in a step-index fiber. 

For the biaxial graded-index fiber the HEn mode is noticeably displaced from the HEn 
modes for the isotropic and uniaxial graded-index fibers. A comparison with figure 3 shows 
that it lies approximately half the distance between the curves for the isotropic step-index 
fiber and the isotropic graded-index fibers. 
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3* Stratification Technique 
3.1 Formulation of Problem 

A purely numerical approach to the problem of solving eq. (2-17) is to subdivide the core 
into N homogeneous layers as shown in figure 14 and then solve an easier problem in each 
layer [17]. Note, this solution method is valid for all modes of an isotropic or uniaxial fiber 
and the transverse modes, i.e. m = 0, for a biaxial graded-index-fiber. For the case m/0 
in a biaxial fiber eqs. (2-17a) and (2-17b) remain coupled and this solution method does 
not work. 


{N+l) 


Figure 14 Geometry for stratification technique 



For the case of an isotropic or uniaxial fiber fi(r) = 62 (f). In the n’th layer eq. (2-17) 
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can be written as 


dr 1 r dr j 

d 2 h[ n) 1 dh[ n) 


(n) 

(n) _ C 3 2 
f 3 (n) K 


nr 


>E { z n) = 0 


dr 2 


+ \ d ^r + {a 2 [4" } - * 2 ] - ^ J = o 


where i — 1, 3, is the approximate value of €{(r) in the n’th layer. If we let 


(n) 

A n ) _ f 3 t .2 

3 >) * 

e i 


= A 

W*)* 


and 


then the solution of eq. (3-1) in the n’th layer is given by 
E { z n) = 

hW = 


A n J m {p[ n)r ) + C n Y m {p[ n) r), 


)’ > 0 

A n I m {q[ n) r) + CnKmiq^r), 

{p ( r y 

) 2 <0 

5 n J m (p ( 2 n) r) + DnYmip^r), 


)’>«* 

B n Im(q[ n) r) + D n K m (q { 2 n) r ), 

«•> 

)’<«* 


(3 - la) 
(3 - 16) 


(3 - 2a) 
(3 - 26) 
(3 - 3a) 

(3 - 36) 


(3 - 4a) 


(3 - 46) 


In order for the fields to be finite in the first layer Ci and D\ must be identically zero. In 
the cladding, jV + 1 layer, we require that Apt+\ = -Bat+i = 0 and 

= (p^ +1 >) 2 = A 2 (t c - «>) < 0 (3 - 5) 

so that the fields are exponentially decaying. 


wr = 


» 


- 


J n ) 

e 3 k 2 

>) K 




Recall that 
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But 





Therefore 

= = (3-6) 

f l 

where /(") = 4^/4”^ Since 4"^ and 4"^ are both positive in the core of the fiber (p^) 2 
and {p[ n) ) 2 will always have the same sign. For (p ( 2 n) ) 2 > 0 the tangential components of 
the electric and magnetic fields are given by 


E ( " ] - A n J m {<J-f^p { 2 n) r) + C n Y Tn (\[f^p { 2 n) r) 
h[ n) = BMp^r) + D n Y m (p[ n) r) 


(3 - 7a) 
(3 - lb) 


rE^ = 


a(4 n) ) 2 

r Bnj;(p (») r) + D„Y^?r)\ 


+ 


a{k\ n) Y 


(3 - 7c) 


(n) _ rn0 

a (4 n> ) 


rh'? 1 = 


BMp^r) + DMp^r) 


~ [AM\ff^P { 2 n) r) + CnY^y/fWpPr)] (3 - Id) 

a {^t ) 2 

where E i"* and h[ n) follow from eq. (3-1), rE^ and rh ^ n) follow from eq. (2-8) and {k[ n) f = 
^ 0 (^ 1 ") _ * 2 )- Note t ^ iat ec l s - ( 3 ’ 7c ) “d (3-7d) can be simplified by using the following 




,( n ) 


N 


3_ _ ,/,<">>> 
(n) - V f l f 3 


(3-8) 


and 





ak[ n ^ 1 

«(*j" ) ) 2 = *F 


(3-9) 
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For (p ( 2 n) ) 2 < 0 the tangential fields are given by eq. (3-7) provided we make the following 


substitutions 


3 m 

Ym - 


4 n) 


K„ 


<£ ] 


V V 

J m 1 m 


Y' 

x m 


K'm 


mfi jkpr jkpr 

a(fc' ( "’) J “('I' 1 "*) 2 ’ Jfe["7 7^ 


10 


where (7* n *) 2 = — (Jbj ) 2 . Eq. (3-7) can be written in matrix form as 


E[ n] \ 

h i n) 

= M n 

-4 n \ 
An 

c„ 



\D n J 


where M n is the chain matrix for the n’th layer and is given by 


M„ = 


d(r) 

0 

*ici(r) 


0 

e i( r ) 

k 2 e 2 (r) 


d i (r) 

0 

fcidi(r) 


-*2\/4 n) 4 n)c 2( r ) k l e l( r ) - fc 2 \Zf ( i n) 4 n)d 2 (»-) 


0 \ 
/l( r ) 

hf2(r) 

*l/l( r )7 


and 


for (p ( 2 n) ) 2 > 0 
ci(r) = J m (v/FWM 
di(r) = Yni^/fWp^r) 
ei(r) = 

h(r) = Ym{p 2 ] r) 

*(r) = 4.(VWM 

d 2 (r) = Y^s/W^r) 
e 2 {r) = JL(P2 n)r ) 

/ 2 (r) = Y^(p[ n) r) 

L _ 

1 a(fc ( , n) ) 2 

**-& 


for (p ( 2 n) ) 2 < 0 
ci(r) = 

di(r) = 

ei(r) = / m (4 n)r ) 

/i( r ) = K m {q^\) 

C2(r) = I' Tn {\/W'<lt ] T) . 
d 2 (r) = K'ni^JW^r) 
e 2 (r) = 7m(4” )r ) 
/2(r) = K^iq^r) 

L _ _ th/9 
Kl _ o(-r(")) J 

*2 = -*& 


(3-10) 


(3-12) 


(3-13) 
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Matching the tangential fields at the surfaces of each layer gives 





/a 2 \ 

M:(rj) 

Bi 

0 

= M 2 (n) 

Bi 

C2 


0 ) 


\Di) 


(A 2 \ 


(A 3 \ 

M 2 (r 2 ) 

B 2 

Cl 

= M 3 (r 2 ) 

b 3 

\Dj 


1 

\D S J 

( 

An \ 


( 


(3-14) 


M N (r N ) 


Bn 
Cn 
\D n ) 


= M N ^(r N ) 


0 

CW+i 

\Dn + i / 


where r n is the normalized radius of the n’th layer. Eq. (3-14) is essentially a system of 4A r 
equations in 4 N unknowns. The propagating modes can be found directly from eq. (3-14) 
by setting the determinant of the system matrix equal to zero. The time required to find a 
determinant of a n x n matrix is proportional to n 3 . If we double the number of layers then 
it takes 8 times as much time to find the determinant. What is needed is a more efficient 
algorithm for determining the allowable modes from eq. (3-14). 

If we recognize that the i’th coefficient vector can be written in terms of the i + l’th 
coefficient vector as 

(t\ (*»\ 

(3 - 15) 

\Z‘J 

then eq. (3-14) can be more conveniently written as 


(Ai\ 


(A l+ i\ 

B t 

Ct J 

= M- 1 (ri)Mi + i(r,-) 

B i+i 

Ci + 1 

\Di) 


Di + 1 / 


Mi(r a ) 


(M\ 
B i 
0 

Vo/ 


= M 2 (r 1 )M- 1 (r 2 )M 3 (r 2 )M 3 - 1 (r3) 


MAr(r;v-i )M N X (r^r )Mjv+i (r^) 


/ 0 \ 

0 

Cn + i 

V Dn + i / 


(3-16) 
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i=2 


Defining an overall chain matrix product, M, where 

N 

M = M 2 (rj) 

eq. (3-16) can then be written as 
hybrid 


(3-17) 


/ 

* \ 







Bx 






Cn+ i 






\£jv+ 1 / 






/(Mi),, 

(Mi),: 

(A4)i3 

~{M) „ \ 

Ai 


(Mi)„ 

(Mi):, 

~(M) 23 

~{MU 

B ! 


(Mi) 31 

(Mi),: 

~{M) 33 

-{MU 

C N+i 


\(Mi)« 

(Mi)« 

-{MU 

-{MU/ 

\Dn + i 


/ 0 \ 
0 
0 

Vo/ 


(3-18) 


where (Mi)^ is an element of M, and {M)^ is an element of M. The problem has been 
effectively reduced from a system of 4 N equations in 4A unknowns into a system of 4 
equations in 4 unknowns. Propagating modes are found by requiring the determinant of 
Mtnd to be identically equal to zero. For the special case of an isotropic or uniaxial step- 
index fiber N = 1 and eq. (3-18) reduces to either eq. (2-29) for the isotropic case or 
eq. (2-31) for the uniaxial case. 

From eq. (3-17) it can be seen that calculating the overall chain matrix, M , requires N- 
1 matrix inversions and 2{N -1) matrix multiplications. Since the size of individual matrices 
is fixed the time needed to find M and hence find the determinant of grows linearly 

with increasing N. The chain matrix approach is therefore a more desirable algorithm for 


determining what modes propagate. 
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3.2 Numerical Results 

The accuracy of this solution method increases with the number of layers, however, with 
this increase in accuracy comes an increase in computation time. For a parabolic profile, a 
choice of five layers gives reasonable accuracy with a minimum of computation time [22] 

Figures 15 and 16 are plots of k versus ko& for the cases m = 0 and 1 of an isotropic 
graded-index fiber with a parabolic permittivity profile. As in the previous example the 
values of n r and n c are taken to be 1.515 and 1.5 respectively. Comparing figures 15 and 
16 with figures 2 and 3 several differences can be seen in the dispersion curves for the step- 
index and graded-index fibers. The most important difference is the value of k as a function 
of ko a increases less rapidly for the graded-index fiber than for the step-index fiber. This 
indicates that the pulse delay which is proportional to dn/d(koa ) is smaller for an isotropic 
fiber with a parabolic permittivity profile than one with a step profile. The other notable 
features are the increase in the cutoff frequencies of all modes except the HEn compared 
with the step-index fiber and several of the hybrid modes have become degenerate or nearly 
degenerate. 

Figures 17 and 18 are plots of k versus k^a for a uniaxial graded-index fiber with 
parabolic permittivity profiles where = 1.515, 713 = 2.0 and n c = 1.5. Comparing figures 
15 and 16 with figures 3 and 4 it can be seen that like the isotropic graded-index fiber the 
value of k. versus k 0 a increases less rapidly than in a uniaxial step-index fiber. The uniaxial 
graded-index also exhibits an increase in the cutoff frequencies for all modes except the 
HEn mode as compared with the uniaxial step-index fiber. 
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4. Discussion 

A total of five types of fiber cores have been discussed. For an isotropic or uniaxial step- 
index fiber either the wave equation formulation or the matrix equation formulation can 
be solved exactly in terms of Bessel functions. The wave equation formulation can be 
solved using WKB analysis for all types except a biaxial graded-index fiber. The matrix 
equation formulation can be solved for all five types of fibers using the method of asymptotic 
partitioning of systems of equations. 

For an isotropic or uniaxial graded-index fiber the WKB solutions are solution of an 
associated scalar wave equation not the full vector problem as given in eqs. (2-17) and (2-18). 
For an isotropic graded-index fiber when the permittivity profile is parabolic and the fields 
are ass um ed to be far from cutoff an approximate solution of the vector problem is possible. 
A comparison of the WKB solutions and the vector solutions shows the vector solutions can 
be obtained from the WKB analysis by renumbering the WKB modes. Strictly speaking this 
comparison is valid only for an isotropic parabolic-index fiber. However, it seem reasonable 
to extend the comparison to isotropic graded-index fibers where the permittivity profiles are 
not parabolic and also to urn axial graded-index fibers. This renumbering of the WKB modes 
has been used as the basis for a more general vector analysis of an isotropic graded-index 
fiber using a generalized WKB technique [10], [15]. 

The negative aspect to the WKB solutions lies with the assumption that the core is 
infinite in extent and therefore there is no need to impose boundary conditions on the 
electric and magnetic fields. For a WKB solution the allowable values of k are determined 
in the process of determining the solution. In contrast, for a step-index fiber, where an 
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exact solution is possible, the allowable values of k are determined by imposing boundary 
conditions. If the renumbered WKB modes are compared with results obtained using the 
stratification technique one finds a poor agreement between the two methods. For example 
according to eq. (2-61) an HEi„ mode is equivalent to a WKBo, n -i mode. However, a 
comparison of figures 6 and 16 shows the HE 12 and HE 13 modes do not correspond to the 
WKB 01 and WKB 0 i modes. The HE i2 and HE i3 modes do agree very well with the WKB 02 
and WKB 04 respectively. This suggests the boundary conditions are important even when 
a mode is far from cutoff. This suggests that further investigation is needed to determine 
whether the WKB modes can be renumbered in such a way as to be valid for a graded-index 
fiber with a finite core and cladding. 

In theory the matrix equation can be solved for any type of fiber core using the method 
of partitioning of systems of equations. The solution obtained is valid wherever the Taylor 
series expansion for the matrix A (s) is valid. Simply because it is theoretically possible to 
solve eq. (2-65) does not mean the solutions obtained are of practical interest. It would be 
useful to compare the asymptotic solution of eq. (2-65) with some known solution, preferably 
an exact solution, in order to determine whether a valid solution has been found. The only 
case where am exact solution is known is for a step-index fiber. 

The asymptotic solution for a step-index fiber can be obtained as a special case of the 
asymptotic solution for either the uniaxial parabolic-index fiber, eq. (2-152) or the biaxial 
graded-index fiber, eq. (2-153). Setting A® and A 3 equal to zero in either eq. (2-152) or 
(2-153) the asymptotic solutions for E z and h z for a step-index fiber can be written as 


Ai(«) 


E z = C l s m P u (s)e 


(4 - la) 
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h z = C 2 3 m P 3 2(^)e Al( * ) 


(4 - 16) 


where 


M *) = “^(«1 ~ K 2 )* 2 


Pn{s) = (Cl - K 2 ) 


1 + 


C3 


4m(m+l)«, (e, "'‘ V ] 

= “ “’'I 1 + 4m(m + l) ( ' 1 " '‘'•’I 


(4 - lc) 
(4 - Id) 
(4 — le) 
(4 - Iff) 


4 m(m + 1) 

and m > 0. Comparing eq. (4-1) with the asymptotic solution of Bessel’s equation given 
by eq. (2-114) E z and h z can be written as 


E z = k 2 Nl Ciy{pk N is) 


h z = k 2 Nl C 2 y{kNis) 


(4- 2a) 


(4 - 26) 


where k^ ll — (i - k 2 , p 2 = e 3 /ei and y(i) is the asymptotic solution for the Bessel function 
Jm and is given by 


»(*) = j* m 


1 + 


€ . 


(4 - 2c) 


4m(m + 1)_ 

Since the exact solution for E z and h z are proportional to JmipkNis) and Jm(kNia) respec- 
tively the asymptotic solutions for E z and h z appear correct. 

In general, when m ^ 0 the solutions to the matrix equation using the method of 
asymptotic partioning of systems of equations appear to be of the general form 


»(*) = P( x ) e 9(x) 


(4-3) 


where p(x) is a monotonic function and q(x) is a positive real monotonically increasing 
function. This implies the solutions will not behave in an oscillatory manner. For example, 
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the asymptotic solution to Bessel’s equation accurately reproduces the form of J m (x) in the 
region Q < z < z\ — 6 where z\ is the first zero of J m and 6 is a small positive number which 
depends on the order of the solution. In the region z > z\ — 6 the asymptotic solution falls 
rapidly to zero and can be taken to be identically equal to zero a short distance past z \ . As 
a consequence only the mode with the lowest cutoff frequency for a given value of m will 
be found when the boundary conditions are imposed. 

This is main disadvantage in solving the matrix equation using the method of asymp- 
totic partitioning of systems of equations. However, the solutions obtained are in good 
agreement with exact and numerical results. On the other hand, the WKB analysis of the 
wave equation produces results which can not be directly compared with either exact or 
numerical results. 

A potental area of further research involves the comparison of the asymptotic solutions 
of the matrix equation with the asymptotic forms of well known functions in an attempt 
either to simplify the problem so as to get better asymptotic solutions or to determine the 


form of the exact solution. 
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Appendix: Computer Programs 

This appendix contains three programs, WKB, ASYMP, and STRAT which were used to 
find the dispersion curves in figures 3-13 and 15-18. The program WKB numerically solves 
the mode condition, eq. (2-43) for either isotropic or uniaxial graded-index fibers where 
the permittivity profiles are given by eq. (2-44). The program ASYMP uses a numerical 
implementation of the method of asymptotic partitioning of systems of equations to solve 
eq. (2-65) and then finds the allowable modes using eq. (2-154). The type of fiber core for 
which ASYMP solves eq. (2-65) is determined by the procedure f indAn. A version of f indAn 
is given for the cases of a uniaxial step-index, a uniaxial graded-index fiber with parabolic 
permittivity profiles and a biaxial graded-index fiber where « 2 ('») is a constant and ei(s) and 
€ 3 ( 3 ) have a parabolic profile. The program STRAT implements the chain matrix version of 
the stratification technique described in section 3. The dispersion curves for the step-index 
fibers, figures 2-5, are obtained using STRAT with the number of layers equal to one. 

All three programs make use of some or all of the following IMSL subroutines: 


bsjO Bessel function of the first kind, order zero 

bsjl Bessel function of the first kind, order one 

bsjs Bessel function of the first kind, real order 

bsyO Bessel function of the second kind, order zero 

bsyl Bessel function of the second kind, order one 

bsys Bessel function of the second kind, real order 

bsiO Modified Bessel function of the first kind, order zero 

bsil Modified Bessel function of the first kind, order one 

bsis Modified Bessel function of the first kind, real order 

bskO Modified Bessel function of the second kind, order zero 

bskl Modified Bessel function of the second kind, order one 

bsks Modified Bessel function of the second kind, real order 

If teg Find LU factorization for a complex general matrix 

If deg Find determinant of a complex general matrix from LU factorization 

merer Multiply two complex rectangular matrices 

ccgcg Copy complex general matrix 

zbren Find zero of a real function which changes sign over a given interval 

zreal Find zeroes of a real function 

qdag Integrate real function using adaptive quadrature 


If any of the above subroutine names is preceeded by the letter d, such as dbsjs, then the 
double precision version is being used. In addition ASYMP has it’s own procedures to perform 
addition, subtraction and multiplication for complex numbers and complex matrices and 
also a procedure to find the determinant of a complex matrix. 


1 C***************** ****************************************** ****** *** *** 

2 C 

3 C This program integerates the phase term of a VKB solution for 

4 C a graded index fiber in order to find the dispersion curve 

5 C 

6 C***************** *********************************** ******************* 

7 C 

8 C Global Constants 

9 C 

10 integer numMax, dSize, maxLvl, maxPnt, iso, uniaxl 

11 real pi 

12 C 

13 parameter ( numMax =5 , dSize=100 ) 

14 parameter ( iso = 1 , uniaxl = 2 ) 

15 C 

16 C Input Parameters 

17 C 

18 C type = 1 for isotropic fiber 

19 2 for uniaxial fiber 

20 C n(l) = maximum value of the refractive index of the core 

21 C in the rho direction 

22 C n(2) = maximum value of the refractive index of the core in 

23 C the phi direction 

24 C nc = refractive index of the cladding 

25 C mu = mode order of the solution 

26 C alf(n) = parameters which describes the shape of the 

27 C refractive index profiles 

28 C Kamin = minimum value of Ka 

29 C Kamax = maximum value of Ka 

30 C numKa = number of divisions between Kamax and Kamin 

31 C KppMin = minimum value of normalized propagation constant 

32 C KppMax = maximum value of normalized propagation constant 

33 C numKpp = number of division between KppMin and KppMax 

34 C 

35 real n(3), alf(3), nc , Kamax, Kamin, KppMax, KppMin 

36 integer mu, numKa, numKpp, type 

37 C 

38 C Computed Parameters 

39 C 

40 C e(i) = maximum value of permittivity in the core 

41 C = n(i)**2 

42 C ec = permittiviy of cladding = nc**2 

43 C delKa = increment for Ka = ( Kamax -Kamin ) /numKa 

44 C delKpp = increment for kappa = ( KppMax-KppMin ) /numKpp 

45 C 

46 real e(3), ec, delKa, delKpp 

47 C 

48 C Program Variables 

49 C 

50 C Ka = ko * a = normalized wave number 

51 C kappa = normalized propagation constant = B/ko 

52 C 

53 C i , j ,k,loopKp,loopKa = loop variables 

54 C 

55 real Ka, kappa, Kpp(dsize) , IntDvPi (dSize ,2) , slope, Ksoln, 

56 ft result , delY 

57 integer i, loopKa, loopKp 

58 logical*l flag(2), modes (2), iroots, uroots 

59 C 

60 C Root finding and integration parameters 

61 C 

62 integer irule 

63 parameter ( irule =2 ) 

64 C 

65 real errabs , errrel, errest, rl(2), r2(2), eps , eta 

66 parameter ( errabs=0 . 001 , errrel=0.001 ,eps=l ,0e-07, eta=0.1 ) 

67 C 

68 C 

69 C Phase function declarations 


75 


70 C 

71 

72 

73 

74 C 

75 C. . 

76 C 

77 

78 C 

79 C. . 

80 C 

81 C 

82 C 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 10 

93 

94 

95 

96 

97 C 

98 

99 
100 
101 
102 

103 

104 

105 

106 

107 11 

108 C 

109 

no 

in 

112 c 

113 

114 

115 C 

116 C 

117 C 

118 

119 

120 

121 C 

122 C 

123 C 

124 

125 

126 

127 C 

128 

129 

130 

131 

132 

133 C 

134 C 

135 C 

136 C 

137 

138 


real psil , psi2 
external psil, psi2 

common kappa, Ka, e, ec, all, mu, type 


pi = acos( -1.0 ) 


Read input parameters 

read (53,*) type 
do 10 i = 1, 3 

read (53,*) n(i), alf(i) 
print 101, i, n(i) , i, alf(i) 
e(i) = n(i)**2 
if ( alf(i) .It. 1 ) then 
print 102, i 
stop 
endif 
continue 

if ( n(l) .ne. n(2) ) then 

print *, ’ >>Error : n(l) and n(2) must be equal’ 
stop 
endif 

read (53,*) mu, nc 
print *, 'nc = * , nc 
print * , 'mu = ' , mu 
ec = nc**2 
do 11 i = 1, 3 

if ( n(i) . le. nc ) then 
print 104, i 
8 top 
endif 
continue 

read (53,*) KaMin, KaMax, numKa 

read (53,*) KppMin, KppMax, numKpp 

if ( numKpp .gt. dsize ) numKpp = dsize 

delKa = (KaMax -KaMin) /numKa 
delKpp = (KppMax-KppMin) /numKpp 

Loop through values of Ka 

delKa = ( Kamax-Kamin ) /numKa 
Ka = Kamin 

do 70 loopKa = 1, numKa+1 

Loop through values of B 

kappa - KppMax + delKpp 
flag(l) = .false, 
flag (2) = .false. 

do 40 loopKp = 1, numKpp 
kappa = kappa - delKpp 
Kpp (loopKp) = kappa 
IntDvPi (loopKp , 1 ) =0.0 
IntDvPi( loopKp ,2) = 0.0 

Find turning points r and r for psil and psi2 


modes (1). = .false . 
modes(2) = .false. 
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160 

161 C 

162 C 

163 C 

164 C 
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185 C 

186 C 

187 C 

188 C 

189 

190 

191 

192 45 

193 

194 

195 

196 

197 

198 

199 

200 
201 
202 

203 

204 50 

205 

206 60 


if (mu . eq. 0 ) then 
do 30 i = 1, 2 

if ( n(i) .ge. kappa ) then 

r2(i) = ( (e(i)-kappa**2)/(e(i)-ec) )**(1 . 0/alf (i)) 
rl (i) = -1.0*r2(i) 
modes (i) = .true, 
else 

rl (i) = 0.0 
r2(i) = 0.0 
endif 
continue 

if ( type .eq. iso ) then 
modes(2) = .false, 
endif 

else 

if ( type .eq. iso ) then 

modes(l) = iroots( rl(l), r2(l), 1 ) 
else 

modes(l) = uroots( rl(l), r2(l) ) 
modes(2) = iroots( rl(2), r2(2), 2 ) 
endif 
endif 

Integrate phase terms from r to r 

1 2 

if ( modes (1) ) then 

call qdag( psil, rl(l), r2(l), errabs, errrel, irule , 
ft result, errest ) 

endif 

IntDvPi (loopKp, 1) = result/pi 
flag(i) = .true, 
else 

IntDvPi(loopKp, 1) = 0.0 
endif 

if ( modes (2) ) then 

call qdag( psil, rl(2), r2(2), errabs, errrel, irule, 
ft result, errest ) 

endif 

IntDvPi(loopKp,2) = result/pi 
flag(i) = .true, 
else 

IntDvPi(loopKp,2) = 0.0 
endif 

continue 

Determine values of kappa for which the integral between the 
turning points satisfies the phase condition. 

do 60 i = 1 , 2 

if ( flag(i) ) then 
delY =0.5 

if ( ( intDvPi(l,i) - delY ) .gt. 0.0 ) then 
delY = delY +1.0 
goto 45 
endif 

do 50 loopKp = 2, numKpp 

if ( ( intDvPi( loopKp, i) - delY ) .gt. 0.0 ) then 

slope = ( intDvPi ( loopKp-1 , i ) - intDvPi ( loopKp ,i) ) 
ft / ( Kpp(loopKp-l) - Kpp(loopKp) ) 

Ksoln = Kpp(loopKp)+(delY-intDvPi (loopKp, i) )/ slope 
delY = delY +1.0 
print 100, i, Ka, Ksoln 
endif 
continue 
endif 
continue 
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207 C 

208 Ka = Ka + delKa 

209 70 continue 

210 C 

211 C 

212 C 

213 100 format( IX, ’>>Phase condition satisfied for psi’, II, 

214 ft ’ at Ka = F6.2, ’ B/Ko = ’, F16.10 ) 

215 101 format( IX, ’n(’, II, ’) = ’, f8.6, ’, alf(’, II, ’) = ’, F5.3 ) 

216 102 formate IX, ’»Error: alf(’, II, ’) must be greater than 0’ ) 

217 103 format (IX, ’»Error: n(’, II, ’) must be greater than nc» ) 

218 C 

219 

220 C 

221 stop 

222 end 

223 C 

225 C 

226 logical function iroots( rl , r2, select ) 

227 C 

228 C****** ************ ******************** ****** *************************** 

990 r 

230 real rl, r2 

231 integer select 

232 C 

233 real delR, delta, errabs , errrel, dOne 

234 integer i, maxfn, num, iso, uniaxl, biaxl 

235 parameter ( delR=0.1, delta=l . Oe-6 , num= 10 ) 

236 parameter ( errabs=0.0, errrel=l . Oe-5 , d0ne=1.0 ) 

237 parameter ( iBo = 1, uniaxl = 2, biaxl = 3 ) 

238 C 

239 real r(num+l), er(num+l), a, b, newSgn, oldSgn 

240 integer count 

241 C 

242 real e(3), alf(3), ec , kappa, Ka 

243 integer mu, type 

244 common kappa, Ka, e, ec, alf , mu, type 

245 C 

246 real gl , g2 

247 external gl , g2 

248 C 

249 

250 C 

251 C Calculate values of function between 0 and 1 

252 C 

253 do 10 i = 1, num+1 

254 r(i) = (i-l)*delR 

255 if ( i .eq. 1 ) then 

256 r(i) = delta 

257 endif 

258 if ( select .eq. 1 ) then 

259 er (i) = gl( r(i) ) 

260 else 

261 er (i) = g2( r(i) ) 

262 endif 

263 10 continue 

264 C 

265 C Look for sign change and then use library routine to find root 

266 C 

267 count = 0 

268 oldSgn = sign( dOne, er(l) ) 

269 do 20 i = 2, num+1 

270 newSgn = sign( dOne, er(i) ) 

271 if ( newSgn .ne. oldSgn ) then 

272 count = count + 1 

273 a = r(i-l) 

274 b = r(i) 

275 maxfn = 15 
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• 297 
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299 
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326 

327 

328 

329 

330 

331 

332 

333 

334 

335 

336 

337 

338 

339 

340 

341 

342 

343 

344 


if ( select .eq. 1 ) then 

call zbren( gl , errabs, errrel, a, b, maxfn ) 
else 

call zbren( g2 , errabs , errrel, a, b, maxfn ) 
endif 

if ( count .eq. 1 ) then 
rl = b ' 
else 
r2 = b 
endif 
endif 

oldSgn = newSgn 
20 continue 
C 

if ( count .eq. 2 ) then 
iroots = .true, 
else 

rl = 0.0 
r2 = 0.0 
iroots = .false, 
endif 
return 
end 
C 

C*+****** **************************************************** *********** 

c 

logical function uroots( rl, r2 ) 

C 

C************************** ********************************************* 

c 

real rl, r2 
C 

real errabs, errrel, dOne 

integer maxfn, num, iso, uniaxl, biaxl 

parameter ( num=10 ) 

parameter ( errabs=0.0, errrel=l . Oe-5 , d0ne=1.0 ) 
parameter ( iso = 1, uniaxl = 2, biaxl = 3 ) 

C 

real r(num+l), er(num+l) , Rmax , delR, delta, a, b, newSgn, oldSgn 
integer count , i 
C 

real e(3), alf(3), ec, kappa, Ka 
integer mu, type 

common kappa, Ka, e, ec, alf , mu, type 
C 

real gl 
external gl 
C 

C 

c 

C Calculate values of function between 0 and 1 
C 

delR =0.1 

if ( type .eq. biaxl ) then 

Rmax = ( (e(2)-kappa**2)/(e(2)-ec) )**(! .0/alf (2)) 
delR = (Rmax-1 . 0e-12)/num 
endif 

delta = 1.0e-5 * delR 
do 10 i = 1, num+1 
r(i) = (i-l)*delR 
if ( i .eq. 1 ) then 
r(i) = delta 
endif 

er(i) = gl( r(i) ) 

10 continue 
C 

C Look for sign change and then use library routine to find root 
C 

count = 0 
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361 
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364 

365 
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367 
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369 
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371 
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373 
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394 C 
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oldSgn = sign( dOne, er(l) ) 
do 20 i = 2, num+1 

nevSgn = sign( dOne, er(i) ) 
if ( newSgn .ne. oldSgn ) then 
count = count + 1 
a = r(i-l) 
b = r(i) 

maxfn =15 ^ . 

call zbren( gl , errabs, errrel, a, b, maxfn ) 


if ( count 
if ( count . 
endif 

oldSgn = newSgn 
continue 


eq. 

eq. 


rl = 
r2 = 


if ( count .eq. 2 ) then 
uroot 8 = .true, 
else 

if ( (count. eq.l) . N x ^ 

fc . and. (type . eq.biaxl) .and. (alf (2) .gt .alf (1) ) ) then 

uroot s = .true. 
r2 = Rmax 
else 

rl = 0.0 
r2 = 0.0 
uroots = .false, 
endif 
endif 
return 
end 
C 

C** ************** ********** 

c 

function psil( r ) 

c**************** ******************************************************* 

c 


real r 

integer iso, uniaxl, biaxl 
parameter ( iso = 1, uniaxl = 

real e(3) , alf (3), ec , kappa, 

integer mu, type 

common kappa, Ka, e, ec, alf, 


2, biaxl = 3 ) 


Ka 


real erl , er2, er3. 


mu, type 
KppSqr , term2 


KppSqr = kappa* *2 

if ( type .eq. iso ) then 

x = e (1 ) - ( e(l) - ec )*abs (r )**alf (1 ) - KppSqr 
else 

erl = e(l) - ( e(l) - ec )*abB(r)**alf (1) 

er2 = e(2) - ( e(2) - ec )*abs (r)**alf (2) 

er3 = e(3) - ( e(3) - ec )*abs(r)**alf (3) 

x = ex3 - er3*(KppSqr)/erl 
endif 

if ( mu .ne. 0 ) then 

term2 = (mu/(Ka*abs(r) ))**2 
if ( type .eq. biaxl ) then 

term2 = term2*er2*(erl-KppSqr)/(erl*(er2-KppSqr) ) 
endif 

x = x - term2 
endif 

if ( x -ge. 0.0 ) then 
psil - Ka*sqrt( x ) 
else 
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psil = 0.0 
endif 
return 
end 


function psi2( r ) 


integer iso, uniaxl, biaxl 

parameter ( iso = 1, uniaxl = 2, biaxl = 3 ) 

real e(3), alf(3), kappa, Ka 
integer mu, type 

common kappa, Ka, e, ec, alf , mu, type 
real er2, x 


414 

415 

416 

417 

418 

419 

420 

421 

422 C 

423 C**************************************** ******************************* 

424 C 

425 real r 

426 C 

427 

428 

429 C 

430 

431 

432 

433 C 

434 

435 C 

436 C. 

437 C 

438 

439 

440 

441 

442 

443 

444 

445 

446 

447 

448 

449 

450 C 

451 c********* ***************************************************** ********* 

452 C 

453 function gl( r ) 

454 C 

455 C************************** ************************************** ******* 

456 C 

457 real r 

458 C 

459 

460 

461 C 

462 

463 

464 

465 C 

466 

467 C 

468 C. 

469 C 

470 

471 

472 

473 

474 

475 

476 

477 

478 

479 

480 

481 

482 


er2 = e(2) - ( e(2) 
x = er2 - kappa**2 
if ( mu .ne. 0 ) then 

x = x - (mu/(Ka*r))**2 
endif 

if ( x .ge. 0.0 ) then 
psi2 = Ka*sqrt( x ) 
else 

psi2 = 0.0 
endif 
return 
end 


ec )*abs(r)**alf (2) 


integer iso, uniaxl, biaxl 

parameter ( iso - 1, uniaxl - 2, biaxl = 3 ) 

real e(3), alf (3), ec , kappa, Ka 
integer mu, type 

common kappa, Ka, e, ec, alf, mu, type 
real erl , er2, er3 , KppSqr, terml, term2 


erl = e (1 ) - ( e(l) - ec )*abs(r)**alf (1) 

er2 = e (2) - ( e(l) - ec )*abs (r)**alf (2) 

er3 = e (3) - ( e(3) - ec )*abs (r)**alf (3) 

KppSqr = kappa* *2 

terml = erl - KppSqr 

term2 = ( mu / (Ka*abs(r)) )**2 

if (type ,ne. iso ) terml = er3*terml/erl 

if ( type .eq. biaxl ) then 

term2 = term2*er2* (erl-KppSqr)/(erl*(er2-KppSqr) ) 
endif 

gl = terml - term2 

return 

end 


483 C 

484 C+++++++***+**+****++**+***+++*+***+***++++******+*+******************** 

485 C 

486 function g2( r ) 

487 C 

488 Q*********************************************************************** 

489 C 

490 real r 

491 C 

492 integer iso, uniaxl, biaxl 

493 parameter ( iso = 1, uniaxl = 2, biaxl = 3 ) 

494 C 

495 real e(3), alf(3), ec , kappa, Ka 

496 integer mu, type 

497 common kappa, Ka, e, ec, alf , mu, type 

498 C 

499 C 

500 C 

501 g2 = e (2) - ( e (2) -ec )*abs(r)**alf (2) - kappa**2 

502 k “ (mu/(Ka*abs(r) ) )**2 

503 return 

504 end 
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1 program Asymp( output ); 

2 

3 const 

4 maxSolnOrder = 10; 

5 maxSize = 100; 

6 mMaxPlus2 = 6; 

7 

8 type 

9 integer2 = -32768 . .32767 ; 

10 complex = record 

11 x, y : real; 

12 end; 

13 modeType = ( HE, EH ); 

14 matrix = array [1 . .4 , 1 . .4] of complex; 

15 superMatrix = array [0 . .maxSolnOrder] of matrix; 

16 eigenvalues = array[1..4] of real; 

17 vector = array [1 . .mMaxPlus2] of real;. 

18 

19 var 

20 noRootsFound : boolean; 

21 i, j, 1, n, m, solnOrder, numKa, numKappa, oldSign, 

22 newSign, debug, colNum : integer2; 

23 nl, el, n3 , e3 , nc , ec , Ka, KaMin, KaMax, kappaMin, kappaMax, 

24 delKappa, delKa, lowerKappa, upperKappa, lowerDet, 

25 upperDet, det , root, oldRoot : real; 

26 cmplxDe term inant : complex; 

27 kappa, determinant : array [1 . .maxSize] of real; 

28 data : text; 

29’ 

30 procedure cmplx( a, b : real; var z : complex ); 

31 begin { cmplx > 

32 z.x := a; 

33 z.y := b; 

34 end; { cmplx } 

35 

36 procedure cmplxAdd( zl, z2 : complex; var result : complex ); 

37 begin { cmplxAdd } 

38 result. x := zl.x + z2.x; 

39 result. y zl.y + z2.y; 

40 end; { cmplxAdd } 

41 

42 procedure cmplxSub( zl , z2 : complex; var result : complex ); 

43 begin { cmplxSub } 

44 result. x := zl.x - z2.x; 

45 result. y := zl.y - z2.y; 

46 end; { cmplxSub } 

47 

48 procedure cmplxMult( zl , z2 : complex; var result : complex ); 

49 var 

50 temp : complex; 

51 begin { cmplxMult } 

52 temp.x := zl.x*z2.x - zl.y*z2.y; 

53 temp.y := zl.x*z2.y + z2.x*zl.y; 

54 result := temp; 

55 end; { cmplxMult > 

56 

57 procedure scalarMult( a: real; z : complex; var result : complex ); 

58 begin { scalarMult } 

59 result. x := a*z.x; 

60 result. y := a*z.y; 

61 end; { scalarMult > 

62 

63 function RealPart( z : complex ) : real; 

64 begin { RealPart } 

65 RealPart := z.x; 

66 end; { RealPart } 

67 
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117 

118 
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123 

124 
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134 


procedure mat Add ( A, B : matrix; var C : matrix ); 
var 

i, j : integer2 ; 
begin { matAdd > 
for i := 1 to 4 do 
for j := 1 to 4 do 

cmplxAdd( A[i,j], B[i,j], C[i,j] ); 
end; { matAdd } 

procedure matSub( A, B : matrix; var C : matrix ); 
var 

i, j : integer2 ; 
begin { matSub } 
for i := 1 to 4 do 
for j := 1 to 4 do 

cmplxSub( A[i,j], B[i,j], C[i,j] ); 
end; { matSub } 

procedure matMult( A, B : matrix; var C : matrix ); 
var 

i, j, k : integer2 ; 
sum, product ; complex; 
begin { matMult > 
for i := 1 to 4 do 
for j := 1 to 4 do 
begin 

cmplx ( 0.0, 0.0, sum ) ; 
for k := 1 to 4 do 
begin 

cmplx Mult ( A [i ,k] , B[k,j], product ); 
cmplxAdd( sum, product, sum ); 
end; 

C [i , j ] := sum; 
end; 

end; { matMult } 

procedure cmplxDet( A : matrix; var determinant : complex ); 
type 

mat 2 = array [1. .2,1. .2] of complex; 
mat3 = array [1. .3,1. .3] of complex; 

var 

subDet , sum, product : complex; 
subMat : mat 3; 

i, j, k, linePnt : integer2; 

procedure cmplxDet2( A : mat2; var determinant : complex ); 
var 

product 1, product2 : complex; 
begin { cmplxDet2 } 

cmplxMult( A [1,1], A [2, 2], product 1 ); 
cmplxMult ( A [1,2], a [2,1], product 2 ); 
cmplxSub( product 1 , product2, determinant ); 
end; { cmplxDet2 > 

procedure cmplxDet3( A : mat3; var determinant : complex ); 
var 

subDet , sum , product : complex ; 
subMat : mat 2; 

i , j , k , linePnt : integer2 ; 
begin { cmplxDet3 > 

cmplx ( 0.0, 0.0, sum ); 
for i := 1 to 3 do 
begin 

linePnt := 1; 
for j := 1 to 3 do 
if ( j <> i ) then 
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begin 

for k := 1 to 2 do 

subMat [linePnt ,k] := A[j,K+l]; 
linePnt : = linePnt + 1 ; 
end; 

cmplxDet2 ( subMat , subDet ) ; 
cmplxMult( A [i , 1] , subDet, product ); 
if odd(i) 

then cmplxAdd( sum, product, sum ) 
else cmplxSub( sum, product, suin ); 

end; 

determinant := sum; 
end; { cmplxDet3 } 

begin { cmplxDet } 

cmplx ( 0.0, 0.0, sum ) ; 
for i := 1 to 4 do 
begin 

linePnt : = 1 ; 
for j := 1 to 4 do 
if ( j <> i ) then 
begin 

for k := 1 to 3 do 

subMat [linePnt ,k] := A[j,K+l]; 
linePnt := linePnt + 1; 
end; 

cmplxDet3( subMat, subDet ); 
cmplxMult( A[i,l], subDet, product ); 
if odd(i) 

then cmplxAdd( sum, product, sum ) 
else cmplxSub( sum, product, sum ); 

end; 

determinant := sum; 
end; { cmplxDet } 

procedure IdentMat ( var A : matrix ); 
var 

cmplxZero, cmplxOne : complex; 
i, j : integer2; 
begin { IdentMat } 

cmplx ( 0.0, 0.0, cmplxZero ); 
cmplx ( 1.0, 0.0, cmplxOne ); 
for i := 1 to 4 do 
begin 

for j := 1 to 4 do 
A[i,j] := cmplxZero; 

A[i,i] := cmplxOne; 
end; 

end; { IdentMat > 

procedure ZeroMat( var A : matrix ); 
var 

cmplxZero : complex ; 
i, j : integer2; 
begin { ZeroMat } 

cmplx ( 0.0, 0.0, cmplxZero ) ; 
for i := 1 to 4 do 
for j := 1 to 4 do 
A [i , j ] := cmplxZero; 
end ; { ZeroMat } 

function power ( x : real; n : integer ) : real; 
var 

i : integer2; 
product : real; 
begin { power } 
product x; 
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for i := 1 to n-1 do 
product := x * product; 
power := product; 
end ; { power > 

{ > 

procedure findPO (m: integer; er, kappa, Ka: real; var P0: matrix); 
var 

f acl , fac2, fac3: real; 
begin { findPO > 


213 

ZeroMat (P0) ; 


214 

fact : = m * kappa; 

215 

fac2 := m * er; 

216 

fac3 := er - kappa * kappa; 

217 

P0[1, 1] .x := 


218 

f ac3; 

219 

P0[1, 3] .x := 

f ac3; 

220 



221 

P0 [2 , l].x := 

f acl ; 

222 

P0[2, 2] .y := 

m; 

223 

P0[2, 3] .x := 

f acl ; 

224 

P0[2, 4] .y : = 

-1.0 * m; 

225 



226 

P0 [3, 2] .x := 

f ac3; 

227 

P0 [3 , 4] .x : = 

f ac3 ; 

228 

P0 [4, 1] .y := 


229 

-f ac2 ; 

230 

P0 [4 , 2] .x : = 

f acl ; 

231 

P0 [4 , 3] .y : = 

f ac2; 

232 

P0 [4 , 4] .x := 

f ac 1 ; 

233 

end; { findPO ] 


234 



235 

procedure findPOinv (m: integer; 

236 


er, kappa, Ka: real; 

237 


var POinv: matrix ); 

238 

var 


239 

facO, f acl , fac2, fac3, fac4, facS: real; 

240 



241 

begin { findPOinv > 

242 

ZeroMat (POinv) ; 

243 

f acl := 1.0 / 

(2 * (er - kappa * kappa)) 

244 

fac2 := kappa 

* facl; 

245 

fac3 := f ac2 / er; 

246 

fac4 := 1.0 / 

(2 * m) ; 

247 

facS := fac4 / er ; 

248 

POinv [1, l].x 


249 

: = facl; 

250 

POinv [1, 3] .y 

: = -1.0 * fac3; 

251 

POinv [1 , 4] .y 

:= facS; 

252 



253 

POinv [2, 1] .y 

:= fac2; 

254 

POinv [2, 2] .y 

: = -1.0 * f ac4; 

255 

POinv [2, 3].x 

:= facl; 

256 

POinv [3 , 1] .x 


257 

:= facl; 

258 

POinv [3, 3] .y 

: = fac3; 

259 

POinv [3, 4] .y 

:= -1.0 * fac5; 

260 

POinv [4, 1] .y 


261 

:= -1.0 * f ac2; 

262 

POinv [4, 2] .y 

: = fac4; 

263 

POinv [4, 3].x 

:= facl; 

264 

end; { findPOinv > 
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{ The implementation of findAn depends upon the type of fiber } 
{ for which the problem is being solved. Separate version of } 
{ findAn for a uniaxial step- index fiber, a uniaxial } 
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{ parabolic-index fiber and a biaxial graded- index fiber can } 
{ be found following this program listing. } 

procedure findAn( n, m : integer 2; 

el, e3, ec, kappa, Ka : real; 
var An : matrix ) ; 


begin { findAn > 
end ; { findAn > 

procedure findBn( Fn : matrix; var Bn : matrix ); 
var 

i : integer2; 
begin { findBn } 

ZeroMat( Bn ); 
for i := 1 to 4 do 
Bn [i , i] := Fn[i,i] ; 
end; { findBn } 


procedure findWn( Fn : matrix; n : integer2; 

evals : eigenvalues; typeOfMode : modeType; 
var Wn : matrix ) ; 


var 

i, j, jLow, jHigh : integer2; 
cmplxZero : complex; 
begin { findWn } 

ZeroMat( Wn ); 
jLow := 1; 
jHigh := 2; 

if ( typeOfMode = EH ) then 
• begin 

jLow := 3; 
jHigh := 4; 
end; 

for i := 1 to 4 do 

for j := jLow to jHigh do 

if ( i <> j ) 

then scalarMult( -1 . 0/ (evals [i] -evals [j] -n) , 

Wn[i, j] ) ; 

end; { findWn } 


Fn[i , j] , 


procedure findFn( POinv : matrix; 

var A, B, P : superMatrix; 
n : integer 2 ; var Fn : matrix ); 

var 

terml , productl, product2, sum : matrix; 

1 : integer2; 
begin { findFn } 

matMult ( A [n] , P [0] , terml ); 
if ( n > 1 ) then 
begin 

ZeroMat( sum ); 
for 1 := 1 to n-1 do 
begin 

matMult ( A [n-1], P[l], productl ); 
matMult ( P [13 , B [n-1] , product2 ); 
mat Add ( sum, productl, sum ); 
matSub( sum, product2, sum ); 
end; 

end; 

mat Add ( terml, sum, sum ); 
matMult ( POinv, sum, Fn ); 
end; { findFn } 


function dbsk0( var x 
fortran; 


real ) : real; 
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function dbskl( var x : real ) : real; 
fortran; 

procedure dbsks( var xnu, x : real; 

var n ; integer2; 
var bsk : vector ) ; 

fortran; 

procedure vsFort; fortran; 

procedure WriteMat( A : matrix ); 
var 

i, j : integer2; 
begin { VriteMat > 
for i := 1 to 4 do 
begin 

for j := 1 to 4 do 

write ( '(>, A[i,j] ,x, A[i,j].y, ’) * ); 

writeln; 
end; 

end; { WriteMat } 

procedure findDet( solnOrder, m : integer2; 

el, e2, e3, kappa, Ka : real; 
var determinant : real ) ; 

var 

typeOfMode : modeType; 

i, n, r, c, col, cLow, cHigh, delC, size, absM : integer2; 
powerOfS, KaSqr, gmmNSq, gmmN, x, Km, KmPrime, xnu : real; 
cmplxDeterminant : complex; 

Vi, bsk : vector; 
evals : eigenvalues ; 
temp : complex; 

A, B, F, P, W : superMatrix; 

POinv, zero, approxP, bndCon : matrix; 

begin { findDet } 
if ( m > 0 ) 

then typeOfMode := HE 
else typeOfMode := EH; 


ZeroMat( zero ); 
for i := 0 to solnOrder do 
begin 


ACi] 

= zero; 

B[i] 

= zero; 

F [i] 

= zero; 

P [i] 

= zero; 

V[i] 

= zero; 

end; 


evalB [1] 

= m; 

evals [2] 

= m; 

evals [3] 

= -l*m; 

evals [4] 

= -l*m; 

for i := 1 to 4 do 


cmplx( evals [l], 


0 . 0 , 


B[0,i,i] 


); 


findP0( m, el, kappa, Ka, P[0] ); 
findP0inv( m, el, kappa, Ka, POinv ); 
findAn( 0, m, el, e3, ec, kappa, Ka, A[0] ); 


if ( solnOrder > 1 ) then 
for n := 2 to solnOrder do 
begin 

findin( n, m, el, e3, ec, kappa, Ka, A[n] ); 
findFn( POinv, A, B, P, n, F[n] ); 
findBn( F[n] , B[n] ); 

findWn( F[n] , n, evals, typeOfMode, V[n] ); 
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matMult ( P[0], V[n] , P[n] ) ; 
end; 

{ Find boundary condition matrix } 

powerOfS := 1; 

KaSqr := Ka * Ka; 
approxP := zpro; 
cLow := 1; 
cHigh := 2; 
delC := 0; 

if ( typeOfMode = EH ) then 
begin 

cLow := 3; 
cHigh := 4; 
delC := “2; 
end; 

for n := 0 to solnOrder do 
if not odd(n) then 
begin 

for r := 1 to 4 do 

for c := cLow to cHigh do 
begin 

scalarMult ( powerOfS, P[n,r,c], temp); 

CmplxAdd( approxP [r , c] , 

temp , approxP [r , c] ) ; 

end; 

if ( n = 0 ) then 
begin 

Vi Cl] := power( Ka, abs(m) ); 

Vi [2] := Vi Cl] ; 
end 
else 
begin 

Vi [1] := Vi[l] * exp( RealPart( B [n,cLow,cLow] ) 

* powerOfS / n ) ; 

Vi [2] Vi [2] * exp( RealPart( B [n, cHigh, cHigh] ) 

* powerOfS / n ) ; 

end; 

powerOfS := KaSqr * powerOfS; 
end; 

bndCon := zero; 
for r := 1 to 4 do 

for c := cLow to cHigh do 
begin 

col := c + delC; 

scalarMult ( Vi[col], approxP [r ,c] , bndCon [r , col] ); 
if not odd( r ) then 

scalarMult ( 1.0/Ka, bndCon[r,c], bndCon[r ,col] ); 

end; 

gmmHSq := kappa * kappa - ec; 
gmmH := sqrt( gmmNSq ); 
x := Ka * gmmN; 
if ( m = 0 ) then 
begin 

Km := dbsk0( x ); 

KmPrime := -1.0 * dbskl( x ); 
end 
else 
begin 

absH := abs( m ); 
size absM + 2; 
xnu := 0.0; 

dbsks( xnu, x, size, bsk ); 

Km := bsk [absM+1] ; 
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KmPrime := -0.6*( b*k[absM] + bsk[absM+2] ); 
end ; . . 

Cmplx( -I .0*Km, 0.0, bndCon[i,3J ); 
bndCon[3,4] := bndCon[l ,3] ; 

Cmplx( m*kappa*Km/ (Ka*gmmHSq) , 0.0, bndConL2,3J 
bndCon[4,4] := bndCon[2,3] ; 

Cmplx ( 0.0, KmPrime / gmmN , bndCon[2,4J ;; 

Cmplx ( 0.0, - 1 . 0 *ec*KmPrime/gmmN, bndCon[4,3j ); 


CmplxDet ( bndCon, cmplxDeterminant ); 
determinant : = cmplxDeterminant .x , 
end; { findDet } 


function sign( x : real ) : integer 2 ; 
begin { sign } 
if ( x <0.0 ) 
then sign := -1 
elBe sign := 1; 
end; { sign > 

function findRoot( xl . yl. x 2 , y 2 : real ) : real 
var 

slope : real; 
begin { findRoot } 

slope : = ( y 2 -yl )/( x 2 ~xl ); 
findRoot : = xl - yl/slope; 
end; { findRoot } 


begin { Asymp } 
vsFort ; 

reset ( data ); 

readln( data, m, solnOrder ); 

reading data, nl , n3, nc ); 

readlnC data, KaMin, KaMax, numKa ) ; 

readlnC data, kappaMin, kappaMax , numKappa ;; 

readlnC data, debug, colNum ); 


el := nl*nl ; 
e3 := n3*n3 ; 
ec := nc*nc; 

delKa := ( KaMax-KaMin ) /numKa; 

delKappa := ( kappaMax -kappaMin ) /numKappa; 

Ka := KaMin; 
for i:= 1 to numKa do 
begin 

for j := 1 to numKappa do 
begin 

kappa [j] := kappaMax - ( j-1) *delKappa; 
findDet ( solnOrder, m, el, el, e3, 

kappa [j] , Ka, determinant LjJ ); 

if ( debug = 0 ) then tin 

writelnC Ka:5:2, kappa [j] : 11 : 6 , determinant [ 3 ] ); 

end; { for j } 

{ Find roots by looking for changes in sign of determinant > 
oldSign := eign( determinant [i] ); 

j • = i » 

noRootBFound := true; . v . 

while ( noRootsFound and ( 3 <= numKappa ) } do 

begin { while } 
j : = j + 1 ; 

newSign := sign( determinant [j] ); 
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if ( newSign <> oldSign ) 
then 
begin 

noRootsFound := false; 
lowerKappa := kappa[j]; 
lowerDet := determinant [j] ; 
upperKappa := kappa [j-1] ; 
upperDet := determinant [j-1] ; 
oldRoot := 1.49; 
root := lowerKappa; 

while ( abs( oldRoot-root ) > 1.0e-6 ) do 
begin 

oldRoot := root; 

root := findRoot( lowerKappa, lowerDet, 

upperKappa, upperDet ) ; 
findDet( solnOrder, m, el, el, e3, 
root , Ka , det ) ; 
if ( det < 0.0 ) 
then 
begin 

lowerKappa := root; 
lowerDet := det; 
end 
else 
begin 

upperKappa := root; 
upperDet := det; 

end; 

end; { while } 

end; 

oldSign := newSign; 
end; { while } 
if noRootsFound 

then writeln( Ka:4:l,* no roots* ) 

else writelnC Ka:4:l, ><’, colNum:l,’> root:10:8); 

Ka := Ka + delKa 
end; { for i } 
end. { Asyrap > 

{ findAn for a step-index fiber } 

procedure findAn ( n, m : integer; 

el, e3, ec, kappa, Ka : real; 
var An : matrix ) ; 


var 


facl, f ac2 , fac3, fac4, fac5, fac6, f ac7 , del 


begin { findAn } 

ZeroMat ( An ) ; 
if ( n - 0 ) then 
begin 

facl := m * kappa; 
fac2 := kappa * kappa; 
fac3 := el - fac2; 
f ac4 := m * m; 


An[l ,3] .y 
An[l ,4] .y 
An [2, 3] .y 
An[2 ,4] .y 
An[3 , 1] .y 
An [3 ,2] .y 
An[4 ,1] .y 
An [4 ,2] .y 


-1 .0*f acl/el 
f ac3/el ; 
f ac4/el ; 
f acl/el ; 
facl ; 

-1 ,0*f ac3 ; 

-1 .0*fac4; 

-1 .0*f acl ; 


end; 


real; 
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if ( n = 2 ) then 
begin 

An [2, 3] .y := -1; 

An [4,1] .y := e3; 
end; 

end ; { findAn } 

{ findAn for an uniaxial graded-index fiber > 

{ where el and e3 have parabolic profiles } 

procedure findAn ( n, m : integer; 

el, e3, ec, kappa, Ka : real; 
var An : matrix ) ; 

var 

facl, fac2, fac3, fac4, fac5, fac6, fac7, del : real; 

begin { findAn } 

ZeroMat( An ); 
if not odd(n) then 
begin 

facl := m * kappa; 
fac2 := kappa * kappa; 
fac3 := el - fac2; 
fac4 := m * m; 
if ( n = 0 ) then 
begin 

An[l , 3] .y := -1 . 0*f acl/el ; 

Anti ,4] .y := fac3/el ; 

An[2,3] .y : = fac4/el; 

An [2 ,4] .y := facl/el ; 

An[3,l] .y := facl; 

An [3 , 2] .y := -1 .0*fac3; 

An[4,l].y := -1.0*fac4; 

An[4,2] .y := -1.0+faci; 
end 
else 
begin 

del :=(el-ec)/(2*el); 
f ac5 := 2 * del / ( Ka * Ka ) ; 
fac6 := power( facS, n div 2 ); 
fac7 := f acl*f ac6/el ; 

An [1 , 3] .y := -1 . 0*f ac7/el ; 

An[2,4] .y := fac7 ; 

An[l,4].y := -1 .0*f ac2*f ac6/el ; 
if (n=2)or(n=4) then 
case n of 
2 : begin 

An[2,3] .y := fac4*fac6/el - 1.0; 
An[3,2].y := el*fac6; 

An [4,1] .y := e3 ; 
end; 

4 : begin 

An[2,3].y := f ac4*f ac6/el ; 

An[4,l].y := ( ec - e3 )/( Ka * Ka ); 
end; 

end 

else An[2,3].y := f ac4*f ac6/el ; 
end; 

end; 

end ; { findAn } 

{ findAn for a biaxial graded-index fiber } 

{ where el and e3 have a parabolic profile } 

{ and e2 is constant } 
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procedure findAn( n, in : integer2; 

el, e3 , «c, kappa, Ka : real; 
var in : matrix ) ; 

var 

facl, fac2, fac3, fac4, fac5, fac6, f ac7 , del : real; 

begin { findAn > 

ZeroMat ( in ) ; 
if not odd(n) then 
begin 

facl := m * kappa; 
fac2 := kappa * kappa; 
fac3 := el - fac2; 
fac4 := m * m; 
if ( n = 0 ) then 
begin 

Anti, 3] .y := -1 .0*f acl/el ; 

An[l ,4] .y := f ac3/el ; 

An [2 ,3] .y := fac4/el ; 

An [2, 4] .y := facl/el ; 

An[3 , 1] ,y := facl ; 

An [3 ,2] .y := -1.0*fac3; 

An [4 , 1] .y := -1 .0*fac4; 

An [4,2] .y := -1 *0*facl ; 
end 
else 
begin 

del := ( el - ec ) / ( 2 * el ) ; 
f ac5 := 2 * del / ( Ka * Ka ) ; 
fac6 := power ( fac5, n div 2 ); 
fac7 := f acl*f ac6/el ; 

An[l, 3] .y := -1 .0*fac7/el ; 

An [2, 4] .y := fac7 ; 

An[l,4].y := -1 . 0*f ac2*f ac6/el ; 
if ( n = 2 ) or ( n = 4 ) then 
case n of 
2 : begin 

An[2,3].y := fac4*fac6/el - 1.0; 

An [4 , 1] .y := e3; 
end; 

4 : begin 

An [2, 3] .y := f ac4*f ac6/el ; 

An[4,l] ,y := ( ec - e3 )/( Ka * Ka ) ; 
end; 

end 

else An[2 ,3] . y := f ac4^f ac6/el ; 
end; 

end; 

end ; { findAn } 


C-SL 


1 c********* ******************************* ******************************* 

2 C 

3 C This program uses an approximate analytical method to calculate 

4 C the dispersion curves for an uniaxial graded index fiber. 

5 C 

6 C********* ******** ******************************* ******** *************** 

Global Constants 

numMax = maximum number of layers 
dSize = size of determinant array 

integer numMax, dSize 

parameter ( numMax=10, dSize=100 ) 

Input Parameters 

nl = maximum value of the refractive index of the core 
in the rho and phi directions 

n3 = maximum value of the refractive index in the z direction 
nc = refractive index of the cladding 
a = radius of the core 
mu = mode order of the solution 

alfl, alf3 = parameters which describes the shape of the 
refractive index profiles 
num = number of layers ( num .le. numMax ) 

KaMin = minimum value of Ka 
KaMax = maximum value of Ka 

numKa = number of divisions between KaMax and KaMin 
KppMax = maximum value of kappa 
KppMin = minimum value of kappa 

numKpp = number of divisions between KppMax and KppMin 

real*8 nl , n3 , sc, alfl, alf3, KaMax, KaMin, KppMax, KppMin, 

6 delKa, delKpp 

integer mu, num, numKa, numKpp 

Computed Parameters 

el 
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28 
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e3 = 


ec = 
dell 
del3 


= maximum value of permittivity of the core in the rho and 
phi directions 
= nl**2 

maximum value of permit ivity of the core in the z 
direction 
= n3**2 

permittiviy of cladding = nc**2 
= ( el-ec )/( 2*el ) 

= ( e3-ec )/( 2*e3 ) 

delKa = increment for Ka = ( KaMax -KaMin ) /numKa 
delKpp = increment for kappa = ( KppMax-KppMin ) /numKpp 
rStep = increment for radius of layers = 1.0/num 
rm = array containing radius for each layer 

eml = array containing permittivity in the rho and phi direction 
for each layer 

cm3 = array containing permittivity in the z direction 

real*8 el, e3, ec, dell, del3, Kastep, rstep, rm(numMax+l) , 
fe eml (numMax+1) , em3(numMax+l) 

Program Variables 

Ka = ko * a = normalized wave number 

kappa = normalized propagation constant in the longitudinal 
direction 

i , j ,k,loopKp,loopKa = loop variables 
real*8 Ka, kappa, K(dsize), D(dSize), nrMin, 
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K oldSgn, newSgn, oldDel, newDel 

integer i, loopKa, loopKp, loopB 

Parameters for finding roots 

real*8 errabs, errrel, a, b, eps , eta, xguess(2), x(2) 

integer maxfn, nroots, itmax , infer (2) 

parameter ( errabs - 0.0, errrel - le-6, maxfn = 10) 

parameter ( eps = 0.0, eta = le-5 ) 

parameter ( nroots = 2, itmax = 10 ) 

Declarations needed to make findD a function of one variable 
so that it can be used with dzbren and dzreal . 

real*8 findD 
external findD 

common eml , em3, rm, Ka, mu, num 


Read input parameters 

read (7,*) num 

read (7,*) nl , alfl 

read (7,*) n3, alf3 

read (7,*) nc, mu 

read (7,*) KaMin,KaMax ,numKa 

read (7,*) KppMin, KppMax, numKpp 

if ( num .gt. numMax ) num = numMax 
if ( numKpp .gt. dsize ) numKpp = dsize 

print *, 'nl = ' , nl , 'profile parameter = alfl = alfl 

print *, 'n3 = n3, 'profile parameter = alf3 = alf3 

print * , 'nc = ' , nc 

print *, 'number of layers = ', num 

print * , 'mu = ' , mu 

Calculate radius of each layer 

rstep = 1.0/num 
rm(l) = rstep 
if ( num .gt. 1 ) then 
do 10 i=2,num 

rm(i)=rm(i-l)+rstep 

continue 

endif 

Calculate values of em 

el = nl**2 
e3 = n3**2 
ec = nc**2 

dell = ( el-ec )/( 2*el ) 

del3 = ( e3-ec )/( 2*e3 ) 

eml(l) = el 

em3(l) = e3 

do 11 i = 2, num 

eml(i) = el*( 1 - 2*dell*rm(i-l)**alf 1 ) 
em3 ( i ) = e3*( 1 - 2*del3^rm(i-l)*^alf 3 ) 
continue 
eml(n\im+l) = ec 
em3(num+l) = ec 

delKa = ( KaMax-KaMin )/numKa 
nr Min = dMinl ( nl , n3 ) 
if ( KppMax .gt. nrMin ) KppMax = nrMin 
delKpp = ( KppMax -KppMin ) /numKpp 


139 C 

140 C 

141 C 

142 

143 

144 C 

145 C 

146 C 

147 

148 

149 

150 

151 

152 40 

153 C 

154 C 

155 C 

156 C 

157 

158 

159 

160 
161 
162 

163 

164 

165 

166 

167 50 

168 C 

169 C 

170 C 

171 

172 

173 

174 

175 

176 

177 

178 

179 


Loop through values of Ka 
Ka = KaMin 

do 70 loopKa = 1, mimKa+1 

Loop through values of B 

kappa = KppMax + delKpp 
do 40 loopKp = 1, numKpp 
kappa = kappa - delKpp 
K( loopKp) = kappa 
D( loopKp) = findD( kappa ) 
continue 

Find roots by looking for a change in the 
sign of the determinant 

oldSgn = dsign( l.OdO, D(l) ) 
do 50 i = 2, numKpp 

newSgn = dsign( l.OdO, D(i).) 
if ( newSgn*oldSgn .It. 0.0 ) then 
a = K (i) 
b = K(i-l) 

call dzbren( findD, errabs, errrel, a, b, maxfn ) 
print 100, Ka, b 
endif 

oldSgn = newSgn 
continue 

Look for closely spaced roots 

oldDel = D (2) - D(l) 
oldSgn = dsign( l.OdO, oldDel ) 
do 60 i = 3, numKpp 
newDel = D(i) - D(i-l) 
newSgn = dsign( l.OdO, newDel ) 
if ( newSgn .ne. oldSgn ) then 

if ( (D (i) *D(i-l) .gt. 0.0) .and. (D(i)*newDel .gt. 0)) then 
xguess(l) = K(i-l) 
xguess(2) = K(i) 


call dzreal( findD, errabs, errrel, eps, eta, nroots, 
fe itmax, xguess, x, infer ) 

print 100, Ka, x(l) 
print 100, Ka, x(2) 
endif 
endif 

oldDel - newDel 
oldSgn = newSgn 
continue 

Ka = Ka + delKa 
continue 


180 
181 
182 

183 

184 

185 

186 

187 

188 60 

189 C 

190 

191 70 

192 C 

193 C. . . 

194 C 

195 100 

196 C 

197 

198 C 

199 stop 

200 end 

201 C 

202 C*********************************************************************** 

203 

204 

205 


format ( IX, , »Root at Ka 


F5.0, J , kappa = } , F10.8 ) 


subroutine ident( mat ) 


206 C********************************************************************* 



207 C 

208 complex*16 mat (4, 4) 

209 integer i, j 

210 C 

211 

212 C 

213 do 92 i = 1, 4 

214 do 91 j = 1, 4 

215 mat ( i , j ) = ( 0.0, 0.0 ) 

216 if ( i.eq.j ) then 

217 mat(i, j) = ( 1.0, 0.0 ) 

218 endif 

219 91 continue 

220 92 continue 

221 return 

222 end 

223 C 

224 C******************************************** ******** ******************* 

225 C 

226 double precision function findD( kappa ) 

227 C 

228 C***************************** ************************* ***************** 

229 C 

230 

231 

232 C 

233 

234 C 

235 

236 

237 

238 C 

239 

240 

241 

242 

243 C 

244 C 


(m) 

Calculate values of transverse wave number p 

2 

KppSq = kappa* *2 
do 15 i = 1, num+1 

pm2(i) = eml(i) - KppSq 
continue 

Find M (r ) 

1 1 

call findM( Ml, rm(l), pm2(l), eml(l), 
k em3(l ) , kappa, Ka, mu, 1, num ) 

-1 -1 

Find product M (r )*M (r )*...*M (r )*M (r ) 

2 12 2 num num-1 num num 

call ident( Mtotal ) 
if ( num . gt . 1 ) then 
do 20 m = 2, num 

Find M (r ) 
m m-1 

call findM( Mm, rm(m~l), pm2(m), eml(m), 
ft em3(m) , kappa, Ka, mu, m, num ) 

call dmcrcr( 4, 4, Mtotal, 4, 4, 4, Mm, 4, 4, 4, prod, 4 ) 
call dccgcgC 4, prod, 4, Mtotal, 4 ) 


245 C 

246 C 

247 C 

248 C 

249 C 

250 

251 

252 

253 15 

254 C 

255 C 

256 C 

257 C 

258 

259 

260 C 

261 C 

262 C 

263 C 

264 C 

265 

266 

267 

268 C 

269 C 

270 C 

271 C 

272 

273 

274 

275 


integer numMax 
parameter ( numMax = 10 ) 

real*8 kappa 

real*8 eml (numMax+1) , em3 (numMax+1 ) , rm(numMax+l ) , Ka 
integer mu, num 

common eml, em3 , rm, Ka, mu, num 

complex*16 Ml(4,4), Mm(4,4), Mtotal(4,4), Mmlnv(4,4), 
& prod(4,4), bndcon(4,4) 

real*8 KppSq, pm2 (numMax+l ) , D 
integer i , j , 1 , m 
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c 

c 

c 

c 

c 


276 

277 

278 

279 

280 
281 
282 

283 

284 

285 20 

286 


Find M 

m 


(r ) 
m 


call findMI( Mmlnv, rm(m) , pm2(m), eml(m), 
t em3(m), kappa, Ka, mu, m, num ) 

call dmcrcr( 4,4, Mtotal, 4, 4, 4, Mmlnv, 4, 4, 4, prod, 4 ) 
call dccgcg( 4, prod, 4, Mtotal, 4 ) 
continue 
endif 

287 

288 

289 

290 

291 

292 

293 

294 

295 

296 

297 

298 

299 

300 

301 

302 

303 30 

304 

305 

306 31 

307 32 

308 

309 

310 

311 

312 C 

313 C**** ********** *********** ********************************************** 

314 C 

315 subroutine findM( M, r, ktNsq, epsl, eps3, kappa, Ka, mu, layer. 


C 

C 

C 

C 


Find M (r ) 
num+1 num 

1 = num+1 

call findM( Mm, rm(num) , pm2(l), eml(l), 
k em3(l), kappa, Ka, mu, 1, num ) 

call dmcrcr( 4, 4, Mtotal, 4, 4, 4, Mm, 4, 4, 4, prod, 4) 
call dccgcg( 4, prod, 4, Mtotal, 4 ) 

Find overall matrix which combines all the boundary 
conditions and find determinant 

do 32 i = 1, 4 
do 30 j =1, 2 

bndcon(i,j) = Ml(i,j) 
cont inue 
do 31 j = 3,4 

bndcon(i,j) = -1 .0*Mtotal(i , j ) 
cont inue 
continue 

call det( bndcon, D ) 

findD = D 

return 

end 


316 £ num ) 

317 C 

318 C**** ************* ************************************************ ****** 

319 C 

320 complex*16 M(4,4) 

321 real*8 r, ktNSq, epsl, eps3, kappa, Ka 

322 integer mu, layer, num 

323 C 

324 real*8 cl, c2, dl , d2, el, e2, fl, f2, x, gmmN, ktN, kl , 

325 k k2, zero, facl, fac2 

326 C 

327 data zero / 0.0 / 

328 C 

329 

330 C 

331 facl = dsqrt( eps3/epsl ) 

332 fac2 = dsqrt( epsl*eps3 ) 

333 if ( mu .eq. 0 ) then 

334 kl = 0.0 

335 else 

336 kl = mu*kappa/(Ka*ktNSq) 

337 endif 

338 if ( ktNSq .gt. 0.0 ) then 

339 ktN = sqrt( ktNSq ) 

340 k2 = r/ ktN 

341 x = Ka*ktN*r 

342 call bessel( facl*x, mu, cl, c2, dl, d2, layer ) 

343 call besseK x, mu, el, e2, fl, f2, layer ) 

344 elBe 


345 

346 

347 

348 

349 

350 

351 

352 

353 

354 

355 

356 

357 

358 

359 

360 

361 

362 

363 

364 

365 

366 

367 

368 

369 

370 

371 

372 

373 

374 

375 

376 

377 

378 

379 


gmmN = sqrt( -1.0*ktNSq ) 
k2 = -1 .0*r / gmmN 
x = Ka*gmmN*r 

call mbesslC facl*x, mu, cl, c2, dl, d2, layer, num ) 
call mbessK x, mu, el, e2, fl, f2, layer, num ) 
endif 

M( 1 , 1 ) = dcmplx ( cl ) 

M (1 , 2 ) = ( 0 . 0 , 0.0 ) 

H(1 ,3) = dcmplx( dl ) 

M ( 1 ,4) = ( 0.0, 0.0 ) 

M(2 ,1 ) = ( 0.0, 0.0 ) 

M(2 ,2) = dcmplx ( el ) 

M(2 ,3) = ( 0.0, 0.0 ) 

M(2 ,4) = dcmplx ( fl ) 

M(3 , 1) = dcmplx ( kl*cl ) 

M(3,2) = dcmplx( zero, k2*e2 ) 

M(3 ,3) = dcmplx ( kl*dl ) 

M(3 ,4) = dcmplx ( zero, k2*f2 ) 

M(4,l) = dcmplx( zero, -1 .0*k2*f ac2*c2 ) 

M(4,2) = dcmplx( kl*el ) 

M(4 ,3) = dcmplx( zero, -1 . 0*k2*f ac2*d2 ) 

M(4,4) = dcmplx( kl*fl ) 

return 

end 


c 

subroutine findMI( M, r, ktNSq, epsl, eps3, kappa, Ka, mu, layer, 
ft num ) 

C 

380 O *♦**♦*♦ ****** ************************************ ****** *************** 

381 C 

382 complex*16 M(4,4) 

383 real*8 r, ktNSq, epsl, eps2, kappa, Ka 

384 integer mu, layer, num 

385 C 

386 

387 C 

388 call f indM( M , r, ktNSq, epsl, eps3, kappa, Ka, mu, layer, num ) 

389 call dlincg( 4, M, 4, M, 4 ) 

390 return 

391 end 

392 C 

394 C 

395 subroutine bessel( x, mu, cl, c2, dl , d2, layer ) 

396 C 

398 C 

real*8 x, cl, c2, dl , d2 
integer mu, layer 


399 

400 

401 C 

402 

403 

404 

405 C 

406 

407 C 

408 C 

409 C 

410 

411 

412 


real* 8 xnu 
integer mum ax 

parameter ( xnu = 0.0, mum ax = 4 ) 
real*8 bs j (mumax+2) , bsy (mumax+2) 


if ( mu . gt . mumax ) then 

print *, , »>Error: mu must be less than or equal to’.mumax 
print *. , »> Program terminated due to error.* 
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413 

414 

415 

416 

417 

418 

419 

420 

421 

422 

423 

424 

425 

426 

427 

428 

429 

430 

431 

432 

433 

434 

435 

436 

437 

438 

439 

440 

441 

442 

443 

444 

445 

446 

447 

448 

449 

450 

451 

452 

453 

454 

455 

456 

457 

458 

459 

460 

461 

462 

463 

464 

465 

466 

467 

468 

469 

470 

471 

472 

473 

474 

475 

476 

477 

478 

479 

480 


stop 

endif 

if ( layer .eq. 1 ) then 
dl = 0.0 
d2 = 0.0 

if ( mu .eq. 0 ) then 
cl = dbsj0( x ) 
c2 = -1.0*dbsjl( x ) 
else 

call dbsj8( xnu, x, mumax+2, bsj ) 
cl = bsj(mu+l) 

c2 = 0.5*( bs j (mu)-bs j (mu+2) ) 
endif 

else 

if ( mu .eq. 0 ) then 
cl = dbsj0( x ) 
c2 = -1.0*dbsjl( x ) 
dl = dbsy0( x ) 
d2 = -1.0*dbsyl( x ) 
else 

call db8js( xnu, x, mumax+2, bsj ) 
call dbsys( xnu, x, mumax+2, bsy ) 
cl = bsj (mu+1) 

c2 = 0 . 5* ( bs j (mu)-bs j (mu+2) ) 
dl = bsy(mu+l) 

d2 = 0 . 5* ( bsy (mu) -bsy (mu+2) ) 
endif 
endif 
return 
end 
C 

O*** **************************************************** *************** 

c 

subroutine mbessl( x, mu, cl, c2, dl , d2, layer, num ) 

C 

C****** ********************************************** ******************* 

c 

real*8 x, cl, c2, dl , d2 
integer mu, layer, num 
C 

real*8 xnu 
integer mumax 

parameter ( xnu = 0.0, mumax = 4 ) 

C 

real*8 bsi(mumax+2) ,bsk (mumax+2) 

C 

C 

c 

if ( mu .gt. mumax ) then 

print ♦, , »>Error: mu must be less than or equal to*, mumax 
print *, , »> Program terminated due to error.* 

stop 
endif 

if ( layer .eq. 1 ) then 
dl = 0.0 
d2 = 0.0 

if ( mu .eq. 0 ) then 
cl = dbsiO( x ) 
c2 = dbsil( x ) 
else 

call dbsis( xnu, x, mumax+2, bsi ) 
cl = b8i(mu+l) 

c2 = 0.5*( bsi(mu)+bsi(mu+2) ) 
endif 
endif 

if (( layer .gt. 1 ) .and. ( layer .le. num )) then 
if ( mu .eq. 0 ) then 
cl = dbsi0( x ) 



481 

482 

483 

484 

485 

486 

487 

488 

489 

490 

491 

492 

493 

494 

495 

496 

497 

498 

499 

500 

501 

502 

503 

504 

505 

506 

507 

508 

509 

510 

511 

512 

513 

514 

515 

516 

517 

518 

519 

520 

521 

522 

523 

524 

525 

526 

527 

528 

529 


c2 = dbsil( x ) 
dl = dbsk0( x ) 
d2 = -1.0*dbskl( x ) 

0^,86 

call dbsis( xnu, x, mumax+2, bsi ) 
call dbsks ( xnu, x, mumax+2, bsk ) 
cl = bsi(mu+l) 

c2 = 0.&*( bsi(mu)+b8i(mu+2) ) 
dl = bBk(mu+l) 

d2 = -0.5*( bsk(mu)+bsk(mu+2) ) 
endif 
endif 

if ( layer .eq. num+1 ) then 
cl = 0.0 
c2 = 0.0 

if ( mu .eq. 0 ) then 
dl = dbsk0( x ) 
d2 = -1.0*dbskl( x ) 
else 

call dbsks ( xnu, x, mumax+2, bsk ) 
dl = bsk(mu+l) 

d2 = -0.5*( bsk(mu)+bsk(mu+2) ) 
endif 
endif 
return 
end 

£****+ ******************* *********************************************** 

c 

subroutine det( mat, D ) 

C* ******************************************* *************************** 

c 

complex*16 mat(4,4) 
real*8 D 
C 

complex*16 f ac (4 ,4) ,detl 
real*8 det2 

integer n,lda,ldaf ac , ipvt (4) 

C V 

parameter ( n = 4, Ida = 4, ldafac = 4 ) 


call dlftcg( n, mat, Ida, fac, ldafac, ipvt ) 

call dlf dcg( n, fac, ldafac, ipvt, detl, det2 ) 

D = dreal( detl * 10.0**det2 ) 

return 

end 
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